Next Article in Journal
Bisimulation for Secure Information Flow Analysis of Multi-Threaded Programs
Previous Article in Journal
A Fast Factorisation of Semi-Primes Using Sum of Squares
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of a Model of Leishmaniasis with Multiple Time Lags in All Populations

School of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623, USA
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2019, 24(2), 63; https://doi.org/10.3390/mca24020063
Submission received: 30 May 2019 / Revised: 13 June 2019 / Accepted: 14 June 2019 / Published: 15 June 2019
(This article belongs to the Section Natural Sciences)

Abstract

:
There are several types of deterministic compartmental models for disease epidemiology such as SIR, SIS, SEIS, SEIR, etc. The exposed population group in, for example SEIS or SEIR, usually represents individuals in the incubation class. Time delays (of which there are several types) when incorporated into a SIR or SIS model, also fulfil the role of the incubation period without necessarily adding another compartment to the model. This paper incorporates time delays into a SIS model that describes the transmission dynamics of cutaneous leishmaniasis. The time lags account for the incubation periods within the sandflies vector, the human hosts and the different animal groups that serve as reservoir hosts. A threshold value, R 0 , of the model is computed and used to study the disease-free equilibrium and endemic equilibrium of the system. Analysis demonstrating local and global stability of the disease-free equilibrium when R 0 < 1 is provided for all n + 1 population groups involved is provided. The existence of an endemic equilibrium is only guaranteed when R 0 > 1 and numerical analysis of the endemic equilibrium for a human host, a vector host and a single animal reservoir host that is globally stable is also provided.
MSC:
Primary: 92B05, 92D25; Secondary: 92D30

1. Introduction

Leishmaniasis is primarily a zoonotic disease caused by protozoan parasites that are transmitted to humans by female sandflies. Animals including dogs, cattle, rodents and many others are the main reservoir of the parasite and occasionally humans in some localities. About 30 different species of sandflies are known to transmit more than 20 species of Leishmania parasites resulting in human infection. The incubation period within humans and animals range from two to eight weeks, while that in the sandflies range from four to 25 days [1]. Several forms of the disease have been identified, with the main types being cutaneous, visceral, and mucocutaneous leishmaniasis. Cutaneous leishmaniasis (CL) is the most dominant form of the disease leading to skin lesions that can lead to lifelong scars and disabilities. Visceral leishmaniasis is the most fatal, causing enlargement of the spleen and liver, while mucocutaneous leishmaniasis affects mostly the mucous membranes of the nose and mouth, destroying them in the process.
There is currently no vaccine for CL and although the disease can be treated, the majority of cases occur in developing countries. Estimates of new cases of CL per year range from about 700 , 000 to 1.2 million or more worldwide [2]. In recent years, CL has become increasingly prevalent in urban areas of Latin America, north Africa and central Asia and the Middle East, with these regions contributing about 95 % of new CL cases annually (see Figure 1). Within the endemic regions, six countries, namely; Afghanistan, Algeria, Brazil, Colombia, Iran and Syria, accounted for more than two thirds of new CL cases in 2015 (see Figure 2) [3]. Symptoms of CL in humans are primarily ulcerated lesions on the skin that can take anywhere between a few months and a few years to heal, though cases lasting longer than a year are rare. The lesions usually leave depressed scars on the skin, and can have debilitating effects depending on where they occur. Treatment when available can hasten the healing of lesions and prevents spreading to other parts of the body. Although recognized as one of the most important and widespread parasitic disease in the world, CL prevention and control remains a challenge for health authorities in some countries [4].
Mathematical modeling and analysis has been at the center of infectious disease epidemiology since the classical works of Ross [5] and Macdonald [6]. Various forms of compartmental models for infectious diseases have been formulated [7]. Several models studying the transmission dynamics of CL have been proposed [8,9,10,11,12,13,14]. See [15,16] for a review of mathematical models for visceral leishmaniasis. The majority of these models are deterministic. Bearing in mind that CL is a zoonotic disease that is transmitted to humans by sandflies, in an attempt to reduce the complexities in analyzing models of the disease that incorporate time delays, some authors have only used a single delay, limiting it either to the human host or animal reservoirs. In [17], a mathematical model for CL that incorporates a time delay between infection and infectiousness for the reservoir (animals) hosts only is provided. Their model does not take into consideration time delays in humans (the incidental) host or sandflies (vector) host. In a more recent paper by Roy et al. [18], a model for CL that focuses on the human host and sandflies vector is provided with a single delay incorporated only in the human host. Because within a single locality many different animals can serve as reservoir for CL, Agyingi et al. [19] developed a susceptible-infectious model that describes the transmission dynamics of cutaneous leishmaniasis. The model incorporated a single vector population, multiple animal populations that serve as reservoirs [20] and a human host population. Because the leishmaniasis parasite undergoes an incubation period within the animal reservoirs, sandflies vector and human host, this paper builds on the work in [19] by incorporating time delays in all populations involved, distinguishing it from previous models. The delays represent the time duration between inoculation of susceptible individuals and them becoming infectious.
The paper is organized as follows: the mathematical model is presented in Section 2, followed by the basic properties for nonnegativity, derivation of a threshold value and the existence of equilibriums in Section 3. The local and global stability of the disease free equilibrium is analytically investigated in Section 4. A numerical study of the endemic equilibrium is considered in Section 5 and concluding remarks are given in Section 6.

2. The Mathematical Model

The model presented below builds on the deterministic model in [19]. The CL model in [19] considers n 1 animal populations that serve as reservoirs, a human host population, and a single sand flies vector population. The model consists of n + 1 susceptible classes and n + 1 infectious classes. In this work, we update all 2 ( n + 1 ) classes of susceptible and infectious populations with time delays that account for the incubation period of the parasite within each population of reservoirs, hosts, and vectors. The schematics of the 2 ( n + 1 ) compartmental SIS model is given in Figure 3. As an example, the susceptible compartment label S 1 interacts with infectious sandflies from the compartment I f , resulting in an outflow rate δ 1 ρ 1 into the compartment I 1 . S 1 is also depleted by natural death at rate σ 1 . The S 1 compartment is populated through a birth rate β 1 N 1 and recovery rate γ 1 I 1 . The same process is repeated in all n + 1 susceptible compartments.
The system of equations associated with the schematic diagram that governs the model are:
d S k d t = β k N k δ k ρ k S k ( t ) N k I f ( t τ k ) + γ k I k σ k S k , k = 1 , , n d I k d t = δ k ρ k S k ( t ) N k I f ( t τ k ) γ k I k σ k I k , k = 1 , , n d S f d t = β f N f k = 1 n δ k ρ f , k I k ( t τ * ) N k S f ( t ) + γ f I f σ f S f d I f d t = k = 1 n δ k ρ f , k I k ( t τ * ) N k S f ( t ) γ f I f σ f I f .
All parameters in the above system of equations are positive and their descriptions are stated in Table 1. All susceptible classes are depleted to populate the infectious classes, and are being replenished at some constant birth ( β ’s) and recovery ( γ ’s) rates. All population groups decay naturally at some constant rate (given by the σ ’s).
We define the initial conditions of the system (1) as
S k ( 0 ) = S k 0 > 0 , I k ( t ) = ϕ k ( t ) > 0 , for k = 1 , , n S f ( 0 ) = S f 0 > 0 , I f ( t ) = φ ( t ) > 0 , t [ τ ¯ , 0 ] , τ ¯ = max { τ 1 , , τ n , τ * } ,
where ϕ k ( t ) for k = 1 , , n and φ ( t ) are continuous functions.
Remark 1.
We note here that assumptions such as cutaneous leishmaniasis not leading to any deaths, and equal birth/death rates have been made in the model (1). In this case, the different total populations, N 1 , , N n , and N f are constant, so that N k = S k ( t ) + I k ( t ) for k = 1 , , n , and N f = S f ( t ) + I f ( t ) . (As a consequence it is clear from (1) that S k ( t ) + I k ( t ) = 0 for k = 1 , , n , and S f ( t ) + I f ( t ) = 0 ).
The following result shows that the model (1) with the initial conditions (2) is well-posed.
Theorem 1.
All solutions of the model (1) with initial conditions (2) are positive and bounded.
Proof. 
Let t * be the smallest positive time value where one of the dependent variables in question turns negative, i.e., Z ( t * ) = 0 and d Z / d t ( t * ) 0 where Z is one of the functions S 1 , , S n , I 1 , , I n , S f , I f . Now this Z can not be any of the S k functions, because d S k / d t ( t * ) > 0 when S k ( t * ) = 0 and t * was the smallest positive time value where any of these functions turn negative, i.e., I k ( t * ) 0 . Similarly, Z can not be any of the I k functions because d I k / d t ( t * ) > 0 as S k ( t * ) > 0 and I f ( t * τ k ) > 0 because of the minimality of t * . Analogous arguments prove that Z can not be S f or I f either, thus the statement is verified.
The boundedness of the solutions follow from the fact that the different populations are constant populations with N k = S k ( t ) + I k ( t ) for k = 1 , , n and N f = S f ( t ) + I f ( t ) . □
Remark 2.
Theorem 1 shows that the system (1) with initial conditions (2) is positively invariant on the region defined by
Ω = { 0 S f , 0 I f , 0 S f + I f = N f , 0 S k , 0 I k , 0 S k + I k = N k , k = 1 , , n } .
Further, by the fundamental theorem of functional differential equations [21], the model admits a unique solution.

3. Analysis of the Model

In this section we derive a threshold value of the model and used it to establish the existence of a positive nonzero equilibrium.
In the analysis that follows, we reduced the 2 ( n + 1 ) system Equations (1) to a smaller system of n + 1 equations by focusing only on the equations for the infectious populations. We achieve this by eliminating the equations for the susceptible populations using the conservation equations N k = S k ( t ) + I k ( t ) for k = 1 , , n and N f = S f ( t ) + I f ( t ) . Note that we have assumed that β k = σ k for k = 1 , , n and β f = σ f . The infective equations of the model are now:
d I k d t = δ k ρ k ( N k I k ) N k I f ( t τ k ) γ k I k σ k I k , k = 1 , , n d I f d t = k = 1 n δ k ρ f , k I k ( t τ * ) N k ( N f I f ) γ f I f σ f I f .
We normalize the above equations by defining the dimensionless variables y = I f / N f and x k = I k / N k for k = 1 , , n . The infective equations now become
d x k d t = c k [ α k ( 1 x k ) y ( t τ k ) x k ] d y d t = c 0 [ ( 1 y ) k = 1 n λ k x k ( t τ * ) y ] ,
where y 0 , x k 0 for k = 1 , , n , c 0 = σ f + γ f and c k = σ k + γ k for k = 1 , , n and where
α k = δ k ρ k N f ( σ k + γ k ) N k , λ k = δ k ρ f , k ( σ f + γ f ) , for k = 1 , , n .
From this point forward, for the purpose of simplicity, we assume that τ 1 = τ 2 = = τ n = τ . It is clear that the normalization leading to the system of the Equations (3) constraints the solution space to the positive unit box B = { ( x 1 , , x n , y ) | 0 x 1 , , x n 1 ; 0 y 1 } . The following result, which is a consequence of Theorem 1, establishes that the solution domain given by the unit box B is positively invariant.
Lemma 1.
All solutions of the system (3) that start in the region B remain in B for all t 0 with the exception of the equilibrium solution at the origin.
The proof of this result follows from Theorem 1.
Next we derive a threshold value of the model (3). It is a useful metric that helps determine whether or not an infectious disease can spread through a population. If the infection is to spread so that we have an outbreak, then we need d x k d t > 0 for k = 1 , , n and d y d t > 0 at t = 0 . Thus we need from (3) that,
α k ( 1 x k ) y ( t τ ) x k > 0 , for k = 1 , , n ( 1 y ) k = 1 n λ k x k ( t τ * ) y > 0 .
The first inequality yields
α k ( 1 x k ) y ( t τ ) > x k ,
for k = 1 , , n , and the second inequality gives
( 1 y ) k = 1 n λ k x k ( t τ * ) > y .
Multiplying both sides of (4) by λ k we get
α k λ k ( 1 x k ) y ( t τ ) > λ k x k ,
for k = 1 , , n , which upon summing leads to
k = 1 n α k λ k ( 1 x k ) y ( t τ ) > k = 1 n λ k x k .
By observing that 1 x k 1 and 1 y 1 , the inequalities in (5) and (6), respectively, become
k = 1 n λ k x k ( t τ * ) > y ,
and
k = 1 n α k λ k y ( t τ ) > k = 1 n λ k x k .
From (7) and (8) we get that
k = 1 n α k λ k y ( t τ ) > y .
Noting again that y ( t τ ) y ( t ) at the start of the infection, the above inequality becomes
k = 1 n α k λ k > 1 .
The quantity on the left hand side gives a threshold value of the system (3), and is denoted by R 0 . Thus, we have
R 0 = k = 1 n α k λ k .
Remark 3.
We remark here that the threshold value R 0 computed above is equivalent to the basic reproduction number obtained in [19]. The basic reproduction number is the average number of secondary cases caused by introducing an infectious individual in a completely susceptible population. Based on the above computation of the threshold value, we see that each animal/human population x k , contributes α k λ k infections towards R 0 .
Further we note that for each population x k , α k is given as the product of the average time spent in an infectious state 1 / ( σ k + γ k ) , the sandflies average biting rate δ k , the transmission rate from human/animal to sandflies ρ k , and the ratio of total sandflies to human/animal population N f / N k . Similarly, λ k is given as the product of the average time spent by sandflies in an infectious state 1 / ( σ f + γ f ) , their average biting rate δ k , and the transmission rate from sandflies to human/animal ρ f , k .
We now turn our attention to computing the equilibriums of the system (3). At the equilibrium points, ( x ¯ 1 , , x ¯ n , y ¯ ) , we have that x ˙ k = 0 , ( k = 1 , , n ) and y ˙ = 0 , leading to the equations
x ¯ k = α k y ¯ 1 + α k y ¯ , ( k = 1 , , n ) and k = 1 n λ k x ¯ k = y ¯ 1 y ¯ .
Eliminating x ¯ k from the above equations yields the following equation:
k = 1 n α k λ k y ¯ 1 + α k y ¯ y ¯ 1 y ¯ = 0 .
The Equation (9) governs the equilibriums of the system and an immediate observation is that y ¯ = 0 and consequently ( x ¯ 1 , , x ¯ n , y ¯ ) = ( 0 , , 0 , 0 ) is an equilibrium point. We call ( 0 , , 0 , 0 ) the disease-free equilibrium (DFE) of the model (3). Further, using the Equation (9), the following result establishes the existence of a unique positive equilibrium other than the DFE, which we refer to as the endemic equilibrium.
Theorem 2
([19]). System (3) has a unique equilibrium solution with positive coordinates if the threshold value R 0 = k = 1 n α k λ k > 1 . If R 0 1 , the system has no equilibria with positive coordinates. The origin is an equilibrium in all cases.
The proof of this result is given in Theorem 3.1 in Agyingi et al. [19].
We remark here that the equations defining the equilibrium points are identical to the equations for the ordinary differential equation model in [19], thus the result describing the equilibriums is also identical. We will investigate the stability of the DFE and the endemic equilibrium in Section 4 and Section 5 respectively.

4. Analysis of the Disease-Free Equilibrium

In this section we study the long term behavior of the disease-free equilibrium of the proposed model. We start by computing the characteristic equation of the delay system (3) which is given as
det ( J 0 + e s τ J τ + e s τ * J τ * s I ) = 0 ,
where I is the ( n + 1 ) × ( n + 1 ) identity matrix and where the Jacobian matrices J 0 , J τ , and J τ * are defined respectively as
J 0 = c 1 α 1 y c 1 0 0 0 0 c 2 α 2 y c 2 0 0 0 c n α n y c n 0 0 0 c 0 k = 1 n λ k x k c 0 ,
J τ = 0 0 0 c 1 α 1 ( 1 x 1 ) 0 0 0 c 2 α 2 ( 1 x 2 ) 0 0 0 c n α n ( 1 x n ) 0 0 0 0 and
J τ * = 0 0 0 0 0 0 0 0 0 0 0 0 c 0 λ 1 ( 1 y ) c 0 λ 2 ( 1 y ) c 0 λ n ( 1 y ) 0 .
Bearing in mind that the origin is the DFE of the system, the characteristic equation evaluated at the DFE yields
det c 1 s 0 0 e s τ c 1 α 1 0 c 2 s e s τ c 2 α 2 0 0 0 c n s e s τ c n α n e s τ * c 0 λ 1 e s τ * c 0 λ 2 e s τ * c 0 λ n c 0 s = 0 ,
and on evaluating the above determinant, the characteristic equation at the disease-free equilibrium becomes
( 1 ) n + 1 j = 0 n ( c j + s ) + ( 1 ) n e s ( τ + τ * ) c 0 k = 1 n j = 1 ; j k n ( c j + s ) c k α k λ k = 0 .
Theorem 3.
The disease-free equilibrium of the system (3) is unstable if R 0 > 1 and stable if R 0 < 1 .
Proof. 
We begin by showing that the disease-free equilibrium is unstable if R 0 > 1 . Denoting the function on the left side of (11) by f ( s ) , that is,
f ( s ) = ( 1 ) n + 1 j = 0 n ( c j + s ) + ( 1 ) n e s ( τ + τ * ) c 0 k = 1 n j = 1 ; j k n ( c j + s ) c k α k λ k ,
when s = 0 , we get
f ( 0 ) = ( 1 ) n + 1 j = 0 n c j + ( 1 ) n j = 0 n c j k = 1 n α k λ k .
Observing that k = 1 n α k λ k = R 0 , the preceding equation becomes
f ( 0 ) = ( 1 ) n j = 0 n c j [ R 0 1 ] .
We see from (12) that if n is odd, then f ( s ) as s , and if R 0 > 1 in (13) then f ( 0 ) < 0 . Similarly, if n is even, then f ( s ) as s , and f ( 0 ) > 0 as long as R 0 > 1 . Thus there exists a positive real root of f ( s ) , establishing the instability of the disease-free equilibrium.
Next we show that the disease-free equilibrium is stable if R 0 < 1 . Re-writing Equation (11) and dividing both sides by ( 1 ) n + 1 , we obtain
j = 0 n ( c j + s ) = e s ( τ + τ * ) c 0 k = 1 n j = 1 ; j k n ( c j + s ) c k α k λ k ,
which is equivalent to
c 0 + s = e s ( τ + τ * ) c 0 c 1 α 1 λ 1 c 1 + s + c 2 α 2 λ 2 c 2 + s + + c n α n λ n c n + s .
Suppose that there is a s 0 which satisfies the above equation, then we have that,
c 0 + s = e s ( τ + τ * ) c 0 c 1 α 1 λ 1 c 1 + s + c 2 α 2 λ 2 c 2 + s + + c n α n λ n c n + s < e s ( τ + τ * ) c 0 c 1 α 1 λ 1 c 1 + c 2 α 2 λ 2 c 2 + + c n α n λ n c n = e s ( τ + τ * ) c 0 k = 1 n α k λ k = e s ( τ + τ * ) c 0 R 0 .
We get from the above computation that
s < c 0 R 0 e s ( τ + τ * ) 1 .
If s 0 and R 0 < 1 , then the righthand side of the inequality (15) is negative which is a contradiction. Therefore all real eigenvalues are negative when R 0 < 1 .
Further, we investigate whether there are complex eigenvalues with positive real parts. Suppose that s = ν + i ω , where ν , ω > 0 , then from Equation (14) we get
c 0 + ν + i ω = e ( ν + i ω ) ( τ + τ * ) c 0 c 1 α 1 λ 1 c 1 + ν + i ω + c 2 α 2 λ 2 c 2 + ν + i ω + + c n α n λ n c n + ν + i ω .
Re-writing the above equation we get
1 + ν c 0 + i ω c 0 = e ν ( τ + τ * ) [ cos ( ω ( τ + τ * ) ) i sin ( ω ( τ + τ * ) ) ] × c 1 α 1 λ 1 ( c 1 + ν ) 2 + ω 2 ( c 1 + ν i ω ) + + c n α n λ n ( c n + ν ) 2 + ω 2 ( c n + ν i ω ) .
Equating the real and imaginary parts of the above equation we respectively obtain
1 + ν c 0 = cos ( ω ( τ + τ * ) ) c 1 α 1 λ 1 ( c 1 + ν ) ( c 1 + ν ) 2 + ω 2 + + c n α n λ n ( c n + ν ) ( c n + ν ) 2 + ω 2 e ν ( τ + τ * ) ω sin ( ω ( τ + τ * ) ) c 1 α 1 λ 1 ( c 1 + ν ) 2 + ω 2 + + c n α n λ n ( c n + ν ) 2 + ω 2 e ν ( τ + τ * ) .
and
ω c 0 = ω cos ( ω ( τ + τ * ) ) c 1 α 1 λ 1 ( c 1 + ν ) 2 + ω 2 + + c n α n λ n ( c n + ν ) 2 + ω 2 e ν ( τ + τ * ) sin ( ω ( τ + τ * ) ) c 1 α 1 λ 1 ( c 1 + ν ) ( c 1 + ν ) 2 + ω 2 + + c n α n λ n ( c n + ν ) ( c n + ν ) 2 + ω 2 e ν ( τ + τ * ) .
Adding the squares of both sides of Equations (16) and (17), we get that
1 + ν c 0 2 + ω 2 c 0 2 = k = 1 n c k α k λ k ( c k + ν ) ( c k + ν ) 2 + ω 2 2 + ω k = 1 n c k α k λ k ( c k + ν ) 2 + ω 2 2 e 2 ν ( τ + τ * ) < k = 1 n c k α k λ k ( c k + ν ) ( c k + ν ) 2 + ω 2 2 + ω 2 k = 1 n c k α k λ k ( c k + ν ) 2 + ω 2 2 < k = 1 n c k α k λ k ( c k + ν ) 2 + ω 2 k = 1 n c k α k λ k ( c k + ν ) 2 2 < k = 1 n c k α k λ k c k 2 + ω 2 k = 1 n c k α k λ k c k 2 2 < k = 1 n α k λ k 2 + ω 2 k = 1 n α k λ k c k 2 < R 0 2 + ω 2 c * 2 R 0 2 ,
where c * = min { c 1 , c 2 , , c n } . The above inequality yields
1 + ν c 0 2 + ω 2 c 0 2 < R 0 2 1 + ω 2 c * 2 ,
which cannot be satisfied by any positive ν and ω values when R 0 < 1 , if we impose the sufficient condition that c * c 0 . Therefore there are no complex eigenvalues with positive real parts. Also observe that if ν = 0 , there are no imaginary eigenvalues since no ω > 0 satisfies the above inequality. We conclude that the disease-free equilibrium is locally asymptotically stable. □
Remark 4.
The condition that c * c 0 in the above proof is only a sufficient condition since the disease-free equilibrium does not exhibit any form of instabilities. The disease-free equilibrium remains stable even when c * c 0 as demonstrated in the next result.
We complete our analysis of the disease-free equilibrium by providing the following global stability result, which is stronger than local stability whenever the threshold value R 0 is smaller than 1.
Theorem 4.
The disease-free equilibrium of the system (3) is globally asymptotically stable if R 0 < 1 .
Proof. 
We start with the first equation of the system (3), that is, for k = 1 , , n we have
d x k d t = c k [ α k ( 1 x k ) y ( t τ ) x k ] .
Recalling from Remark 1 that all solutions are contained within a unit box, we have the differential inequalities
d x k d t c k [ α k y ( t τ ) x k ] , k = 1 , , n ,
which are the same as
d x k d t + c k x k c k α k y ( t τ ) k = 1 , , n .
Multiplying both sides of the preceding inequalities by e c k t , we obtain that
d x k d t e c k t + c k x k e c k t = d d t ( x k ( t ) e c k t ) c k α k e c k t y ( t τ ) , k = 1 , , n .
Integrating the above inequalities on the interval ( 0 , t ) , we obtain
x k ( t ) e c k t x k ( 0 ) 0 t c k α k e c k s y ( s τ ) d s , k = 1 , , n .
Rearranging these inequalities, we have that
x k ( t ) e c k t x k ( 0 ) + e c k t 0 t c k α k e c k s y ( s τ ) d s , k = 1 , , n .
Now this implies that for k = 1 , , n ,
x k ( t ) e c k t x k ( 0 ) + e c k t sup y ( t ) 0 t c k α k e c k s d s = e c k t x k ( 0 ) + sup y ( t ) ( α k α k e c k t ) .
Thus we get that
lim sup t x k ( t ) α k lim sup t y ( t ) , k = 1 , , n .
In a similar fashion, taking the second equation of the system (3), that is,
d y d t = c 0 [ ( 1 y ) k = 1 n λ k x k ( t τ * ) y ] ,
andbearing in mind Remark 1, we obtain the differential inequality
d y d t c 0 [ k = 1 n λ k x k ( t τ * ) y ] ,
which is the same as
d y d t + c 0 y c 0 [ k = 1 n λ k x k ( t τ * ) ] .
Again, multiplication by e c 0 t gives
d y d t e c 0 t + c 0 y e c 0 t = d d t ( y ( t ) e c 0 t ) c 0 e c 0 t [ k = 1 n λ k x k ( t τ * ) ] ,
andthen integration on the interval ( 0 , t ) yields
e c 0 t y ( t ) y ( 0 ) 0 t c 0 e c 0 s λ 1 x 1 ( s τ * ) d s + + 0 t c 0 e c 0 s λ n x n ( s τ * ) d s .
We again rewrite this as
y ( t ) e c 0 t y ( 0 ) + e c 0 t 0 t c 0 e c 0 s λ 1 x 1 ( s τ * ) d s + + e c 0 t 0 t c 0 e c 0 s λ n x n ( s τ * ) d s .
Taking the supremums we get
sup y ( t ) e c 0 t y ( 0 ) + e c 0 t sup x 1 ( t ) 0 t c 0 e c 0 s λ 1 d s + + e c 0 t sup x n ( t ) 0 t c 0 e c 0 s λ n d s .
After integration, we then obtain
lim sup t y ( t ) λ 1 lim sup t x 1 ( t ) + + λ n lim sup t x n ( t ) .
Combining the inequalities (18) and (19), we get
lim sup t y ( t ) α 1 λ 1 lim sup t y ( t ) + + α n λ n lim sup t y ( t ) = lim sup t y ( t ) k = 1 n α k λ k .
The preceding inequality is the same as
lim sup t y ( t ) R 0 lim sup t y ( t ) .
Since R 0 < 1 , we get that lim sup t y ( t ) = 0 ; consequently, by (18), lim sup t x k ( t ) = 0 for k = 1 , , n . This concludes the proof. □

5. Numerical Results and Discussion

In this section, we continue our analysis of the equilibriums by considering the case n = 2 of the system (3), where the x 1 variable represents a proportion of infectious human population, x 2 variable represents a proportion of infectious animal population and y is a proportion of infectious sandflies. Here we numerically confirm the theoretical results established in the previous section on the stability of the disease-free equilibrium. Next we investigate the stability of the endemic equilibrium for different parameter values of the model, in particular, the population densities and time delays.
The parameter values used in the analysis are mostly similar to parameter values in related works [17,18,19], or arbitrarily chosen to explore the behavior of the model. The average biting rates of human and animal by sandflies are equal and set at δ 1 = δ 2 = 0.25 per day. The transmission rate from human and animal to sandflies are also equal and set at ρ 1 = ρ 2 = 0.25 . The transmission rate from sandflies to humans and animals are equal and set at ρ f , 1 = ρ f , 2 = 0.25 . The death rate of humans, animals and sandflies respectively are σ 1 = 0.0001 per day, σ 2 = 1 / 365 per day, and σ f = 1 / 14 per day. Values for the recovery rates of humans, animals and sandflies respectively are γ 1 = 12 / 365 per day, γ 2 = 5 / 365 per day and γ f = 1 / 14 per day.
Because animals are mostly the reservoir of the parasite, the initial conditions for all simulations reported below were set at x 1 ( 0 ) = 0 , x 2 ( 0 ) = 0.1 and y ( 0 ) = 0 . We performed multiple simulations in which the initial conditions were perturbed and no changes in the behavior of the equilibrium solutions were observed.
In the simulations that follow, we use two sets of time delays, { τ = 28 , τ * = 14 } which is well within the known incubation periods, and extreme values { τ = 1000 , τ * = 1000 } which are biologically impossible. The very extreme set of time delays was chosen simply to investigate the nature of the stability of the endemic equilibrium, that is, whether there exist critical delay values at which bifurcations take place.
The calculated threshold value R 0 of the system for the stated parameter values is given as R 0 = α 1 λ 1 + α 2 λ 2 , or R 0 = 0.2985 N f / N 1 + 0.5988 N f / N 2 . It is therefore evident that the sandflies to human/animal ratios determines whether R 0 < 1 or R 0 > 1 . We examine each of these cases below.
We begin with the case where the ratios N f / N k 1 for k = 1 , 2 and where R 0 < 1 . The time evolution of all infectious classes is demonstrated in Figure 4, where there are fewer sandflies compare to humans/animals, and in Figure 5, where the densities of sandflies, humans and animals are equal. For population densities N 1 = N 2 = 5000 and N f = 2500 , the results in Figure 4a are for the time delays { τ = 28 , τ * = 14 } , while the results in Figure 4b are for the time delays { τ = τ * = 1000 } . Similarly, for equal population densities N 1 = N 2 = N f = 12 , 500 , the results in Figure 5a are for the time delays { τ = 28 , τ * = 14 } , while the results in Figure 5b are for the time delays { τ = τ * = 1000 } . The results in Figure 4 and Figure 5 show that the disease dies out in both cases since R 0 < 1 . The presence of very large delays only induces small oscillations whose amplitude decreases with time and all solutions converge to the disease-free equilibrium that is asymptotically stable.
We now turn our attention to the case where for k = 1 , 2 , the ratios N f / N k > 1 and where R 0 > 1 . First we note that it is possible for all k = 1 , 2 , the ratios N f / N k > 1 to still lead to R 0 < 1 . As an example, if N f = 5200 , and N 1 = N 2 = 5000 then we get R 0 = 0.9332 < 1 , leading to the stability of the DFE. Since the endemic equilibrium only exists if R 0 > 1 , we focus on those cases where at least one of the ratios N f / N k 1 for k = 1 , 2 and that R 0 > 1 .
The simulations in Figure 6 give the solution profiles of the model for population densities N 1 = N 2 = 5000 and N f = 12 , 500 . The computations in Figure 6a are for the time delays { τ = 28 , τ * = 14 } , while those in Figure 6b are for the time delays { τ = τ * = 1000 } . The results indicate the global stability of the endemic equilibrium. Given very large time delays as illustrated in Figure 6b, the endemic equilibrium at the onset induces some form of instability that grows with time and on attaining maximum amplitude begins to fissile out, and then eventually becomes stable. Although the human and animal population densities are equal, the results in Figure 6 show that the disease is more prevalent in the animal reservoir. This can be attributed to the initial conditions.
The next results given in Figure 7 are for population densities N 1 = 5000 , N 2 = 12 , 500 and N f = 12 , 500 . As in previous computations, Figure 7a is for time delays { τ = 28 , τ * = 14 } and Figure 7b is for time delays { τ = τ * = 1000 } . The simulations in Figure 7 also indicate the global stability of the endemic equilibrium. Exceedingly large time delays that initially rattle the solutions (see Figure 7b), do not render the endemic equilibrium unstable with time. We also observe here that when the animal (reservoir) population is significantly higher than that of human, then the model predicts higher disease prevalence among humans.
Finally, we consider the case with population densities N 1 = 12 , 500 , N 2 = 2500 and N f = 12 , 500 . Note that in this case, the reservoir density is significantly lower than that of human. The simulations given in Figure 8 indicate a very high level of disease prevalence in the reservoir population and the least disease prevalence among humans. This is unlike in all other cases indicated above where the least prevalence has always been among sandflies.

6. Conclusions

In this paper we have presented a model with time delays for the transmission dynamics of cutaneous leishmaniasis, an infectious disease whose prevalence is on the increase. Although mostly located in tropical regions, global warming may provide suitable conditions for its vector to migrate to subtropical regions and spread the disease. The life cycle of the parasite begins in animals which are usually the reservoirs and manifest as some form of the disease in humans via sandflies transmission.
Mathematical models for leishmaniasis pale in comparison to other infectious diseases and as noted before, only a handful of such models account for the incubation period of the parasite within reservoirs, vectors and humans. The novelty of the model constructed here is that time delays serving as the incubation period have been incorporated into all population groups, which has often been avoided or neglected in past works. Starting with an existing deterministic SIS model, delays were inserted into the 2 ( n + 1 ) dimensional system of equations. Because all populations involved were assumed to be constant, the dimension of the model studied was reduced to a system of only n + 1 infective equations by eliminating the susceptible terms.
A threshold value R 0 of the model was computed as a sum of the products of the infected sandflies and the resulting human/animal infections for each human/animal population group. We used R 0 to analyze equilibriums of the model. The disease-free equilibrium of the model is both locally and globally stable when R 0 < 1 . A numerical study of the positive endemic equilibrium for the case n = 2 which only exist when R 0 > 1 , shows that it is globally asymptotically stable even when very extreme time delay values are employed.
From the model simulations we observe that as long as the sandflies population is kept lower or near the level of the human and animal populations, the disease will die out. We also learn that if the sandflies population is substantially higher than that of humans and/or animals, the disease will persist in all populations. In the situation where the disease is persistent, the model predicts the following: (i) the prevalence is higher in animals for equal human and animal populations; (ii) the prevalence is higher in humans when the human population is smaller than that of animals; and (iii) the prevalence is least in humans when the animal population is smaller than the human population.
The main limitation of the model presented in this work is the assumption that all infected hosts will survive the incubation period. A possible way of addressing that is the inclusion of an exponential decay term that represents the average proportion of infected individuals which survive the incubation period. The threshold value R 0 then will depend on the incubation periods as well.

Author Contributions

Conceptualization: E.A and T.W; data curation: E.A and T.W.; methodology: E.A.and T.W; writing—original draft: E.A. and T.W.; writing—review and editing: E.A. and T.W.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Weiss, F.; Vogenthaler, N.; Franco-Paredes, C.; Parker, S.R. Leishmania tropica-induced cutaneous and presumptive concomitant viscerotropic leishmaniasis with prolonged incubation. Arch. Dermatol. 2009, 145, 1023–1026. [Google Scholar] [CrossRef] [PubMed]
  2. Center for Disease Control and Prevention. Available online: https://www.cdc.gov/parasites/leishmaniasis/epi.html (accessed on 14 June 2019).
  3. WHO. Global Health Observatory Data Repository. Available online: http://apps.who.int/gho/data/node.main.NTDLEISHCNUM?lang=en (accessed on 14 June 2019).
  4. Hotez, P.J.; Bottazzi, M.E.; Paredes, C.F.; Ault, S.K.; Periago, M.R. The neglected tropical diseases of Latin America and the Caribbean: A review of disease burden and distribution and a roadmap for control and elimination. PLoS Negl. Trop. Dis. 2008, 2, e300. [Google Scholar] [CrossRef] [PubMed]
  5. Ross, R. The Prevention of Malaria; E. P. Dutton and Company: New York, NY, USA, 1910. [Google Scholar]
  6. Macdonald, G. The measurement of malaria transmission. Proc. R. Soc. Med. 1955, 48, 295–298. [Google Scholar] [PubMed]
  7. Brauer, F. Mathematical epidemiology: Past, present, and future. Infect. Dis. Model. 2017, 2, 113–127. [Google Scholar] [CrossRef] [PubMed]
  8. Barradas, I.; Caja Rivera, R.M. Cutaneous leishmaniasis in Peru using a vector-host model: Backward bifurcation and sensitivity analysis. Math. Methods Appl. Sci. 2018, 41, 1908–1924. [Google Scholar] [CrossRef]
  9. Bacaer, N.; Guernaoui, S. The epidemic threshold of vector-borne diseases with seasonality. The case of cutaneous leishmaniasis in Chichaoua, Morocco. J. Math. Biol. 2006, 53, 421–436. [Google Scholar] [CrossRef] [PubMed]
  10. Chaves, L.F.; Hernandez, M. Mathematical modelling of American Cutaneous Leishmaniasis: Incidental hosts and threshold conditions for infection persistence. Acta Tropica 2004, 92, 245–252. [Google Scholar] [CrossRef] [PubMed]
  11. Lysenko, A.J.; Beljaev, A.E. Quantitative approaches to epidemiology. In The Leishmaniasis in Biology and Medicine; Peters, W., Killick-Kendrick, R., Eds.; Academic Press: London, UK, 1987; Volume I, pp. 263–290. [Google Scholar]
  12. Rabinovich, J.E.; Feliciangeli, M.D. Parameters of Leishmania braziliensis transmission by indoor Lutzomyia ovallesi in Venezuela. Am. J. Trop. Med. Hyg. 2004, 70, 373–382. [Google Scholar] [CrossRef] [PubMed]
  13. Burattini, M.N.; Coutinho, F.A.B.; Lopez, L.F.; Massad, E. Modelling the dynamics of leishmaniasis considering human, animal host and vector populations. J. Biol. Syst. 1998, 6, 337–356. [Google Scholar] [CrossRef]
  14. Davies, C.R.; Llanos-Cuentas, E.A.; Pyke, S.D.M.; Dye, C. Cutaneous leishmaniasis in the Peruvian Andes: An epidemiological study of infection and immunity. Epidemiol. Infect. 1995, 114, 297–318. [Google Scholar] [CrossRef] [PubMed]
  15. DebRoy, S.; Prosper, O.; Mishoe, A.; Mubayi, A. Challenges in modeling complexity of neglected tropical diseases: A review of dynamics of visceral leishmaniasis in resource limited settings. Emerg. Themes Epidemiol. 2017, 14, 10. [Google Scholar] [CrossRef] [PubMed]
  16. Bi, K.; Chen, Y.; Zhao, S.; Kuang, Y.; Wu, C.J. Current Visceral Leishmaniasis Research: A Research Review to Inspire Future Study. BioMed Res. Int. 2018, 2018, 9872095. [Google Scholar] [CrossRef] [PubMed]
  17. Das, P.; Mukherjee, D.; Sarkar, A.K. Effect of Delay on the Model of American Cutaneous Leishmaniasis. J. Biol. Syst. 2007, 15, 139–147. [Google Scholar] [CrossRef]
  18. Roy, P.K.; Biswas, D.; Basir, F.A. Transmission Dynamics of Cutaneous Leishmaniasis: A Delay-Induced Mathematical Study. J. Med. Res. Dev. 2015, 4, 11–23. [Google Scholar]
  19. Agyingi, E.; Ross, D.; Bathena, K. A model of the transmission dynamics of leishmaniasis. J. Biol. Syst. 2011, 19, 237–250. [Google Scholar] [CrossRef]
  20. Ashford, R.W. Leishmaniasis reservoirs and their significance in control. Clin. Dermatol. 1996, 14, 523–532. [Google Scholar] [CrossRef]
  21. Kuang, Y. Delay Differential Equations with Applications in Population Dynamics; Academic Press: San Diego, CA, USA, 1993. [Google Scholar]
Figure 1. Number of new cases of cutaneous leishmaniasis within the most endemic region of the world between 2005–2015. Data from the World Health Organization (WHO) [3].
Figure 1. Number of new cases of cutaneous leishmaniasis within the most endemic region of the world between 2005–2015. Data from the World Health Organization (WHO) [3].
Mca 24 00063 g001
Figure 2. Number of new cases of cutaneous leishmaniasis for the most endemic countries of the world between 2005–2015. Data from WHO [3].
Figure 2. Number of new cases of cutaneous leishmaniasis for the most endemic countries of the world between 2005–2015. Data from WHO [3].
Mca 24 00063 g002
Figure 3. Schematic diagram of the SIS model. The dash lines represent interaction between infectious and susceptible populations, while the solid lines represent population movement into and out of a compartment. For i = 1 , 2 , , n , f , the susceptible compartment S i interacts with infectious sandflies from the compartment I f , with an outflow rate δ i ρ i into compartment I i . S i decreases by a natural death rate of σ i . The S i compartment is populated through births at a rate of β i N i and through recovery at a rate of γ i I i .
Figure 3. Schematic diagram of the SIS model. The dash lines represent interaction between infectious and susceptible populations, while the solid lines represent population movement into and out of a compartment. For i = 1 , 2 , , n , f , the susceptible compartment S i interacts with infectious sandflies from the compartment I f , with an outflow rate δ i ρ i into compartment I i . S i decreases by a natural death rate of σ i . The S i compartment is populated through births at a rate of β i N i and through recovery at a rate of γ i I i .
Mca 24 00063 g003
Figure 4. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 5000 and N f = 2500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Figure 4. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 5000 and N f = 2500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Mca 24 00063 g004
Figure 5. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 12500 and N f = 12500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Figure 5. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 12500 and N f = 12500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Mca 24 00063 g005
Figure 6. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 5000 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Figure 6. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = N 2 = 5000 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Mca 24 00063 g006
Figure 7. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = 5000 , N 2 = 12 , 500 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Figure 7. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = 5000 , N 2 = 12 , 500 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Mca 24 00063 g007
Figure 8. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = 12500 , N 2 = 2500 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Figure 8. Time evolution of infectious human population x 1 ( t ) , infectious animal population x 2 ( t ) and infectious sandflies population y ( t ) for parameter values δ 1 = δ 2 = 0.25 , ρ 1 = ρ 2 = 0.25 , ρ f , 1 = ρ f , 2 = 0.25 , σ 1 = 0.0001 , σ 2 = 1 / 365 , σ f = 1 / 14 , γ 1 = 12 / 365 , γ 2 = 5 / 365 , γ f = 1 / 14 , N 1 = 12500 , N 2 = 2500 and N f = 12 , 500 . The plots given in (a) are for τ = 28 and τ * = 14 , while the plots in (b) are for τ = τ * = 1000 .
Mca 24 00063 g008
Table 1. Definition of the parameters of the model.
Table 1. Definition of the parameters of the model.
ParameterDefinition
S k susceptible human/animal population for k = 1 , , n
I k infectious human/animal population for k = 1 , , n
N k total human/animal population for k = 1 , , n
S f susceptible sandfly population
I f infectious sandfly population
N f total sandfly population
β k birth rate for human/animal for k = 1 , , n
β f sandfly reproduction rate
σ k natural death rate for human/animal for k = 1 , , n
σ f sandfly death rate
γ k recovery rate for human/animal for k = 1 , , n
γ f sandfly recovery rate
δ k average biting rate of human/animal by sandflies for k = 1 , , n
ρ k transmission rate from human/animal to sandfly for k = 1 , , n
ρ f , k transmission rate from sandfly to human/animal for k = 1 , , n
τ k incubation period in human/animal population for k = 1 , , n
τ * incubation period in sandfly population

Share and Cite

MDPI and ACS Style

Agyingi, E.; Wiandt, T. Analysis of a Model of Leishmaniasis with Multiple Time Lags in All Populations. Math. Comput. Appl. 2019, 24, 63. https://doi.org/10.3390/mca24020063

AMA Style

Agyingi E, Wiandt T. Analysis of a Model of Leishmaniasis with Multiple Time Lags in All Populations. Mathematical and Computational Applications. 2019; 24(2):63. https://doi.org/10.3390/mca24020063

Chicago/Turabian Style

Agyingi, Ephraim, and Tamas Wiandt. 2019. "Analysis of a Model of Leishmaniasis with Multiple Time Lags in All Populations" Mathematical and Computational Applications 24, no. 2: 63. https://doi.org/10.3390/mca24020063

APA Style

Agyingi, E., & Wiandt, T. (2019). Analysis of a Model of Leishmaniasis with Multiple Time Lags in All Populations. Mathematical and Computational Applications, 24(2), 63. https://doi.org/10.3390/mca24020063

Article Metrics

Back to TopTop