A Deterministic Model of Schistosomiasis with Spatial ... - CiteSeerX
A Deterministic Model of Schistosomiasis with Spatial Structure Fabio Augusto Milner, Ruijun Zhao Department of Mathematics, Purdue University, West Lafayette, IN 47907-2067 Abstract It has been observed in several settings that schistosomiasis is less prevalent in segments of river with fast current than in those with slow current. Some believe that this can be attributed to flush-away of intermediate host snails. However, free-swimming parasite larvae are very active in searching suitable hosts, which indicates that the flush-away of larvae may also be very important. In this paper, the authors establish a model with spatial structure that characterizes the density change of parasites following the flush-away of larvae. It is shown that the reproductive number, which is an indicator of prevalence of parasitism, is a decreasing function of the river current velocity. Moreover, numerical simulations suggest the mean parasite load is low when the velocity of river current flow is sufficiently high.
1
Introduction
Schistosomiasis, a parasite (Schistosome)-induced disease, is also known as bilharziasis after Theodor Bilharz, who first identified the parasite in Egypt in 1851. Infection is widespread with a relatively low mortality rate, but a high morbidity rate, causing severe debilitating illness. According to WHO [10], an estimated 170 million people in Sub-Saharan Africa, and a further 30 million in North Africa, Asia and South America, suffered from schistosomiasis in 2004. The disease is generally associated with rural poverty. Schistosomes have to go through an intermediate host (snails in most cases) to complete their life-cycle: eggs, miracidia, cercariae, to the final stage of adult flukes. There have been many studies in mathematical modeling using unstructured models [1, 2, 5, 11], and age-structured models [4, 9]. However, none of them incorporated spatial structure in the habitat, such as a river system. Indeed, the river current plays an important role in spreading schistosomiasis, either by affecting intermediate host snails or parasite larvae. It has been noticed in [7] that the probability for explorers in Africa to be infected in the Southern segment of the Omo River is larger than in the Northern part. Two 1
main features in the Southern Omo River are believed to contribute to the increased incidence, slow flow and large number of residences. In [8] Utzinger, et al. have observed that certain snails (such as B. pfeifferi) have preference for a certain range of river current velocities, which may explain the low incidence of schistosomiasis due to lack of snails. Fingerut and his colleagues, [6], found that environmentally triggered downward swimming can quickly bring larvae (Planktonic cercariae) to the bed, promoting contact with benthic intermediate hosts under most flow conditions, which explains the prevalence of infected snails. The main aim in this paper is to investigate how the river current can change the dynamics of schistosomiasis by considering spatial structure only in parasite larvae, namely miracidia and cercariae. The basic structure of the mathematical model is still as in [1, 5, 11], with the following assumptions: (H1) Human beings grow logistically in the absence of parasites; (H2) Snail population grows logistically in the absence of parasites; (H3) Infected Snails do not reproduce; (H4) Distribution of parasites among human hosts is overdispersed (negative binomial); (H5) Birth rates of hosts are greater than their death rates, bj −µj0 > 0, j = h, s. It is natural to incorporate spatial structure for parasite larvae first of all, since miracidia and cercariae are free-swimming in river. However, for convenience and in order to have a simpler model, human hosts, parasites and snails will be considered spatially independent. Since human beings are active more or less everywhere in their habitat, their contribution to the environment can be viewed as approximately homogeneous. The same is true for parasites since they reside inside the human body. The present paper is structured as follows: in the next section, a mathematical model with spatial structure is given; in Section 3, its well-posedness is established; in Section 4, a reduced model—spatially uniform, given by a dynamical system—is studied analytically and numerically; in Section 5, the full spatial model is partly analyzed and simulations are provided. Finally, we present some discussion and draw some conclusions in Section 6.
2
A Mathematical Model with Spatial Structure
Schistosomiasis is usually prevalent in villages near a river or a lake. For simplicity, we assume the aquatic habitat to be a river [0, L], where 0 is the origin of the river and L represents the end of the interested segment of river. Thus, miracidia and cercariae move along with the river current with velocity v. Moreover, we assume that there is no parasitism at the origin of the river, while at the end of the river, the parasitism is determined by the system. 2
In [11] it was shown that the long-term exponential growth of humans used in [5] is not adequate, not only based on demographic considerations, but also on epidemiological ones since the mean parasite load stays the same irrespective of treatment rate unless such rate is sufficient to kill all parasites. In this paper, we focus on the logistic growth in the human population, which is more adequate for long-term demographic projection and it also does away with the problem of mean parasite load not diminishing with treatment. Similar results are expected if we change the logistic model to one with constant recruitment. For simplicity, the velocity of the river flow is assumed to be constant. It has been observed that the miracidia and cercariae move in the river with directional preference—e.g. towards to light [6], and towards intermediate hosts snails. In this paper, we neglect these anisotropic movement and assume that the larvae move only due to the river movement. In other words, the spatial structure is reflected on the assumption that the parasite larvae follow a first order hyperbolic equation with convection coefficient exactly equal to the velocity of river current. Prior to building the model, let us introduce some demographic and epidemiological parameters. bh µh0 µh1 Lh α β λ κ µp r bs µs0 µs1 Ls ρ ds bp ζ µm µc v
= = = = = = = = = = = = = = = = = = = = =
per capita birth rate of human beings per capita natural death rate of human beings additional per capita death rate of human beings due to overcrowding carrying capacity of environment for human beings parasite-induced death rate of human beings per parasite production rate of cercariae per infected snail attachment rate of cercariae to human beings “clumping parameter” of binomial distribution of adult parasites per capita death rate of parasites treatment-induced rate of parasites per capita birth rate of snails per capita natural death rate of snails additional per capita death rate of snails due to overcrowding carrying capacity of environment for snails contact rate (per snail) of uninfected snails with miracidia parasite-induced death rate of snails mean number of eggs (miracidia) laid per parasite “efficacy” of cercariae in infecting humans (0 ≤ ζ ≤ 1) per capita death rate of miracidia per capita death rate of cercariae river current velocity
Let H(t), P (t), S(t), I(t) denote, respectively, the total numbers of human hosts, adult parasites, uninfected snails and infected snails at time t. Let m(x, t), c(x, t) denote, respectively, the density of free-living miracidia and cercariae at Z L Z L time t and location x, and C(t) = c(x, t)dx, M (t) = m(x, t)dx denote, 0
0
3
respectively, the total numbers of miracidia and cercariae at time t. The mathematical model we propose consists of the following system of differential equations:
dH dt
= (bh − µh0 − µh1 H)H − αP,
dP dt
=
dS dt
= (bs − µs0 − µs1 (S + I))S − ρM S,
dI dt ∂m ∂t ∂c ∂t
ζHC − α
¡ κ+1 ¢ P 2 κ
H
− (µh0 + µh1 H + µp + α + r)P,
(1) =
ρM S − (µs0 + µs1 (S + I) + ds )I,
=
−v
∂m + bp P − µm m, ∂x
=
−v
∂c + βI − µc c ∂x
along with homogeneous Dirichlet boundary conditions at the inflow boundary for the two partial differential equations, m(0, t) = 0,
c(0, t) = 0.
(2)
Remark 2.1. The singularity at H = 0 in System (1) can be avoided by introducing two new variables, the mean parasite load per person R = P/H, and total population size of snails T = S + I. It is easy to derive the following equivalent system to System (1) describing the evolution of H, R, S, T, m, c:
dH dt
=
(bh − µh0 − µh1 H)H − αRH,
dR dt
=
ζC −
dS dt
=
(bs − µs0 − µs1 T )S − ρM S,
dT dt ∂m ∂t ∂c ∂t
α 2 R − (bh + µp + α + r)R, κ
(3) =
(bs − µs0 − µs1 T )T − (bs + ds )(T − S)
=
−v
∂m + bp RH − µm m, ∂x
=
−v
∂c + β(T − S) − µc c. ∂x 4
Remark 2.2. Free-living miracidia and cercariae can live up to 48 hours in fresh water, and they are actively searching hosts only in the first several hours. Therefore, we can neglect the small portion of those penetrating into hosts and assume the loss of larvae is only due to natural death. Remark 2.3. Homogeneous Dirichlet boundary conditions in this model are assumed for simplicity and they can be replaced by measured field data when a specific region needs to be considered.
3
Well-Posedness
Concerning the well-posedness of System (1), we can equivalently investigate System (3). Since (3.v) and (3.vi) are first order hyperbolic equations, we can solve them exactly along characteristic lines x = τ + vt, where τ is a parameter. Denote m(t) ˜ = m(τ + vt, t); then Eq. (3.v) is equivalent to dm ˜ + µm m ˜ = bp RH, dt
(4)
and also to
d µm t (e m) ˜ = bp eµm t RH. (5) dt Since the habitat has finite length L, we can group characteristic lines as above x = vt and below x = vt. In the region {(x, t) : vt < x < L}, integrating Eq. (5) from 0 to t, we have Z µm t
e
t
m(t) ˜ = m(0) ˜ + bp
eµm s R(s)H(s)ds,
(6)
0
and then Z
t
m(x, t) = e−µm t m(x − vt, 0) + bp
e−µm (t−s) R(s)H(s)ds.
(7)
0
x Similarly, in the region {(x, t) : 0 < x < vt}, integrating Eq. (5) from t − to v t, we have Z t x x eµm t m(t) ˜ = eµm (t− v ) m(t ˜ − ) + bp eµm s R(s)H(s)ds, (8) x v t− v and then, m(x, t) =
− µvm x
e
Z =
bp 0
x v
x m(0, t − ) + bp v −µm ( x v −s)
e
Z
t
e−µm (t−s) R(s)H(s)ds
t− x v
x x R(s + t − )H(s + t − )ds. v v 5
(9)
In summary, we can write the solution as Z t −µm t e m(x − vt, 0) + b e−µm (t−s) R(s)H(s)ds if x > vt, p 0 m(x, t) = Z xv x x x bp if x ≤ vt. e−µm ( v −s) R(s + t − )H(s + t − )ds v v 0
(10)
Integrating Eq. (10) in the spatial variable x on the interval [0, L], the total Z L population of miracidia M (t) = m(x, t)dx at time t is given by 0
Z vt Z xv x x x e−µm ( v −s) R(s + t − )H(s + t − )dsdx if t < b p v v 0 0 ¸ Z t Z L· + e−µm t m(x − vt, 0) + bp e−µm (t−s) R(s)H(s)ds dx, vt 0 x Z Z L v x x x bp e−µm ( v −s) R(s + t − )H(s + t − )dsdx if t ≥ v v 0 0
L , v (11) L . v
Similarly, we can write the solution of Eq. (3.vi) as Z t −µc t c(x − vt, 0) + β e−µc (t−s) (T (s) − S(s))ds if x > vt, e 0 Z xv c(x, t) = x x x β e−µc ( v −s) (T (s + t − ) − S(s + t − ))ds if x ≤ vt, v v 0 Z L and the total population of cercariae C(t) = c(x, t)dx at time t as
(12)
0
Z L· ¸ Z t −µc t −µc (t−s) e m(x − vt, 0) + β e (T (s) − S(s))ds dx vt 0 Z vt Z xv x x x e−µc ( v −s) (T (s + t − ) − S(s + t − ))dsdx if t < +β v v 0 0 Z L Z xv x x x β e−µc ( v −s) (T (s + t − ) − S(s + t − ))dsdx if t ≥ v v 0 0
L v,
(13)
L v.
Concerning the positivity of solutions to System (3), subtracting Eq. (3.iii) from Eq. (3.iv), we have d(T − S) = (bs − µs0 − µs1 T )(T − S) − (bs + ds )(T − S) + ρM. dt
(14)
System (3) implies that H and S are nonnegative if the initial data are nonnegative. It is clear from Eq. (14) that if M is nonnegative, then T − S 6
is nonnegative. Moreover, from Eq. (3), if T − S is nonnegative, then R is nonnegative. Eq. (11) shows that M is nonnegative if R is nonnegative, and Eq. (13) shows that C is nonnegative if T − S is nonnegative. Therefore, when H(0) > 0 and S(0) > 0, if either R(0) > 0 or T (0) > S(0), then for small ε > 0 we certainly have H, R, S, M, C and T −S positive on (0, ε). It is clear that Eq. (3.i) gives the positivity of H, Eq. (3.iii) gives the positivity of S, and then Eq. (3.iv) gives the positivity of T as long as they exist. Hence, it suffices to show that T − S must remain positive as long as it exists, since that is sufficient (in fact, equivalent) to establishing the positivity of R. Suppose that T − S vanishes for some positive time, and let t¯ be the smallest such time. Then, T − S is positive on [0, t¯) and T (t¯) = S(t¯). It follows that C(t¯) > 0 from Eq. (13). Then, R is positive on [0, t¯] from Eq. (3.ii) and then Eq. (14) implies that T − S is positive on [0, t¯], contradicting T (t¯) = S(t¯). Thus T − S remains positive as long as it exists, and so do R, M and C. Moreover, m(x, t) and c(x, t) are also positive if R(t), H(t), T (t) − S(t) are positive. This completes the proof of positivity of solutions to System (3). Concerning their boundedness, it is evident that Eq. (3.a) and Eq. (3.d) implies that H ≤ Lh and T ≤ Ls . Moreover, S ≤ Ls and I ≤ Ls . For convenience, let us introduce a notation δ = bh + µp + α + r. It is not hard to show that C is bounded by observing that for t > L/v, Z |C(t)| =
L
Z
x v
β 0
x
e−µc ( v −s) (T (s + t −
0
Z
L
≤
βLs µc
=
· ¸ µc βLs v L + (e− v L − 1) µc µc
,
¯ C.
x x ) − S(s + t − ))dsdx v v
x
(1 − e−µc v )dx (15)
0
Notice from Eq. (3.ii) that α dR = ζC − R2 − δR ≤ ζ C¯ − δR, dt κ
(16)
The boundedness of R follows from R ≤ R0 e−δt +
ζ C¯ ζ C¯ ¯ (1 − e−δt ) ≤ R0 + , R. δ δ
7
(17)
Last, for the boundedness of M , it is sufficient to check that for t > L Ã /v, Z L Z xv x x x |M (t)| = bp e−µm ( v −s) R(s + t − )H(s + t − ))dsdx v v 0 0 ≤
¯Z L x bp Lh R (1 − e−µm v )dx µm 0
=
¸ ¯· bp Lh R v − µc L L+ (e v − 1) µm µm
,
¯. M
(18)
Now let us introduce a vector V¯ = (H, R, S, T )t . Then System (3) can be rewritten as dV¯ = F¯ (V¯ ) = (f1 , f2 , f3 , f4 )t , (19) dt where f1 (H, R, S, T ) =
(bh − µh0 − µh1 H)H − αRH,
f2 (H, R, S, T ) =
ζC −
α 2 R − (bh + µp + α + r)R, κ
f3 (H, R, S, T ) = (bs − µs0 − µs1 T )S − ρM S, f4 (H, R, S, T ) =
(bs − µs0 − µs1 T )T − (bs + ds )(T − S),
and M and C are functions of T − S, H, and R defined in (11) and (13). It is easy to show that F¯ is Lipschitz continuous and Picard Iteration Theorem (fixed-point theory) proves local existence and uniqueness of the solution. The global existence and uniqueness of solutions is guaranteed by their a priori boundedness. Then, the following theorem has been proved. Theorem 3.1. System (1), or equivalently, System (3) admits a unique solution for all time. Moreover, if H(0) and S(0) are positive and either T (0) > S(0), or M (0) > 0, or R(0) > 0, or C(0) > 0, then H, R, S, m and c are positive and T > S as long as they exist.
4
Special Case: Model in Spatially Uniform Distribution
At the first glance, it seems strange to talk about a case where all species are uniformly distributed in space since the homogeneous in-flow boundary conditions for miracidia and cercariae will force a disease free solution. However, spatially uniform distribution can happen at steady state if non-homogeneous in-flow boundary conditions are imposed. 8
The diversity of steady states for the uniformly distributed model has been established in [11], except for β in [11] replaced by βL/µc and bp by bp L/µm . In[11] the population size of miracidia is assumed to be simply proportional to that of parasites and the population size of cercariae is proportional to that of infected snails. It is clear that System (3) also admits four parasite-free steady states E0 E1 E2 EDF E
= (H0 , R0 , S0 , T0 , m0 , c0 ) = (H1 , R1 , S1 , T1 , m1 , c1 ) = (H2 , R2 , S2 , T2 , m2 , c2 ) = (HDF E , RDF E , SDF E , TDF E , mDF E , cDF E )
= (0, 0, 0, 0, 0, 0), = (0, 0, Ls , Ls , 0, 0), = (Lh , 0, 0, 0, 0, 0), = (Lh , 0, Ls , Ls , 0, 0).
where, EDF E represents (H, P, S, I, m, c) = (Lh , 0, Ls , 0, 0, 0) for System (1). E0 is the trivial equilibrium, while E1 and E2 are semi-trivial in that they represent, respectively, the ultimate state of a logistic population of snails and of humans in the absence of parasitism. Concerning the non-trivial equilibria, as in [11], we see from System (3) that at steady state H satisfies the cubic equation H 3 + a2 H 2 + a1 H + a0 = 0, where
a2
=
−Lh ,
a1
=
rs α αµs1 µi + ¯ ¯bp )2 , ρ¯bp µh1 κζ β(ρ
a0
=
−
(20)
δαµs1 µi rh αµs1 µi ¯ h1 (ρ¯bp )2 − ζ βµ ¯ h1 (ρ¯bp )2 . κζ βµ
¯ µi are defined as Here, ¯bp , β, ¯bp = bp L , β¯ = βL , µi = bs + ds . µm µc Define the discriminant of the above cubic polynomial as D = Q3 + R2 , where Q=
3a1 − a22 , 9
U=
9a1 a2 − 27a0 − 2a32 . 54
Note that a2 < 0, a1 > 0, and a0 < 0, and we have a similar result as in [11]. Theorem 4.1. Let D, Q and U be defined as above, and consider D and U as functions of treatment rate r. • If Q > 0 or if Q < 0, U (0) > 0 and D(0) > 0, then Eq. (20) admits only one positive root; • If Q < 0, and if D(0) < 0, then Eq. (20) admits three, then two, and finally one positive root as r increases. 9
• If Q < 0, U (0) < 0 and D(0) > 0, then Eq. (20) admits one, two, then three, two, and finally one positive root as r increases. Remark 4.1. It can be easily seen fromµSystem ¶ (3) that if a root H of Eq. (20) r H ρ¯bp rh Lh h satisfies 0 < H < Lh and rs > ρ¯bp 1− H (or rs > ), then α Lh 4α System (3) or, equivalently, System (1) has a positive equilibrium with human population size H.
Bifurcation Diagram 1200
Carrying capacity of human Lh
1000 IV 800 III 600
400
II
200
0
I 0
0.2
0.4 0.6 treatment rate r
0.8
1
Figure 1: Bifurcation Diagram. In region I, there is no interior positive equilibria and EDF E is globally asymptotically stable. In region II, there is one and only one interior positive equilibrium and it is globally asymptotically stable. In region III, there are three positive interior equilibria and two of them are locally asymptotically stable and no oscillatory solution happens. In region IV, there are also three positive interior equilibria and two of them are unstable and oscillatory solution occurs for certain initial data. Figure 1 shows the diversity of system (1) in the number of positive interior equilibria with respect to two bifurcation parameters Lh and r. The parameters are chosen as: L = 1, µh0 = 1/70, bh = 3/140, α = 0.00001, ζ = 0.000027, κ = 0.243, µp = 0.2, µs0 = 0.5, bs = 80, Ls = 20000, ρ = 0.00004, ds = 6.0, µm = 180, µc = 180, β = 4000µc , bp = 20µm . Remark 4.2. Indeed, Figure 1 is a bifurcation diagram of System (1) for the 10
5
Human Hosts
population size
1.2
625
624.5
624
0 6
population size
8
x 10
500 time (in year) Larvae
4 2
0
Parasites
1 0.9
0
500 time (in year) Snails
1000
2000 M C
6
0
x 10
1.1
0.8
1000
population size
population size
625.5
500 time (in year)
1000 500 0
1000
S I
1500
0
500 time (in year)
1000
Figure 2: Oscillatory solution when r = 0 and Lh = 800 on [0, 1000].
5
Human Hosts
population size
1.2
625.45
625.4
625.35 990 6
population size
8
x 10
995 time (in year) Larvae
Parasites
1.1 1 0.9
995 time (in year) Snails
1000
2000 M C
6 4 2 0 990
x 10
0.8 990
1000
population size
population size
625.5
995 time (in year)
1000
S I
1500 1000 500 0 990
995 time (in year)
1000
Figure 3: Oscillatory solution when r = 0 and Lh = 800 on [990, 1000].
11
chosen parameters. We have been unable to derive a stability analysis for positive equilibria. Numerical simulations show that System (1) may converge to a steady state or oscillatory solutions. Figure 2 and Figure 3 show that an oscillatory solution occurs when parameters are in region IV, where the simulations were done using r = 0 and Lh = 800. However, concerning the stability of boundary steady states, we obtain here similar results as in [11], too. Theorem 4.2. E0 , E1 and E2 are unstable. Proof. For convenience, introduce another notation: κ+1 . κ The Jacobian of System (3) at the equilibrium E = (H, R, S, T, m, c) is a11 −αH 0 0 0 0 0 a22 0 0 0 ζL 0 0 a −µ S −ρLS 0 33 s1 . J2 (E) = 0 (bs + ds ) a44 0 0 0 b p R bp H 0 0 −µm 0 0 0 −β β 0 −µc α ¯=α
where
Then,
J2 (E0 ) =
rh 0 0 0 0 0
a11
=
rh − 2µh1 H − αR,
a22
=
−
a33
=
rs − µs1 T − ρLm,
a44
=
rs − 2µs1 T − (bs + ds ).
0 −δ 0 0 0 0
2α R − δ, κ
0 0 rs (bs + ds ) 0 −β
0 0 0 0 0 0 −(µs0 + ds ) 0 0 −µm β 0
0 ζL 0 0 0 −µc
.
It is obvious that λ = rh > 0 is an eigenvalue of J3 (E0 ); then E0 is unstable. Next, rh 0 0 0 0 0 0 −δ 0 0 0 ζL 0 0 0 −µ L −ρLL 0 s1 s s . J2 (E1 ) = 0 (bs + ds ) −rs − (bs + ds ) 0 0 0 0 0 0 0 −µm 0 0 0 −β β 0 −µc 12
Again λ = rh > 0 is an eigenvalue of J3 (E1 ) and then E1 is unstable. Finally, −rh −αLh 0 0 0 0 0 −δ 0 0 0 ζL 0 0 rs 0 0 0 J2 (E2 ) = 0 0 (b + d ) −(µ + d ) 0 0 s s s0 s 0 b p Lh 0 0 −µm 0 0 0 −β β 0 −µc
.
Note that λ = rs > 0 is an eigenvalue of J3 (E2 ) and then E2 is unstable. Regarding the stability of the Disease Free Equilibrium, we have the following result. Define so-called reproductive number of parasitism ¯ ¯ ¯ 0 = ζ βρbp Lh Ls . R (21) δµi ¯ 0 < 1, then EDF E is locally asymptotically stable; if R ¯0 > 1 Theorem 4.3. If R then EDF E is unstable. Proof. Note that System (1) and System (3) are equivalent. It turns out to be easier to find the threshold in this case by looking at the Jacobian of System (1) than System (3). The Jacobian of System (1) at the equilibrium (H, P, S, I, m, c) is rh − 2µh1 H −α 0 0 0 0 a21 a22 0 0 0 ζLH 0 0 a −µ S −ρLS 0 33 s1 , J1 (E) = 0 0 ρLm − µs1 I a44 ρLS 0 0 bp 0 0 −µm 0 0 0 0 β 0 −µc where
µ a21
= ζLc + α ¯ µ
2P H
P H
¶2 − µh1 P,
¶
a22
= −¯ α
− (µh0 + µh1 H + µp + α + r),
a33
= rs − µs1 (2S + I) − ρLm,
a44 = −(µs0 + µs1 (S + 2I) + ds ). At the equilibrium EDEF , the Jacobian is −rh −α 0 0 0 −δ 0 0 0 0 −r −r s s J1 (EDF E ) = 0 0 0 −(b + ds ) s 0 bp 0 0 0 0 0 β 13
0 0 −ρLLs ρLLs −µm 0
0 ζLLh 0 0 0 −µc
.
The characteristic polynomial of the above Jacobian is f (λ) = (λ + (bh − µh ))(λ + (bs − µs )) (22) [(λ + δ)(λ + µi )(λ + µm )(λ + µc ) − bp ρLs ζβLh L2 ]. It is clear that λ = −(bh − µh ) and λ = −(bs − µs ) are two roots of f (λ) and they are negative from our assumptions. Whether f (λ) has positive roots is completely determined by the forth order multiplier f1 (λ)
=
(λ + δ)(λ + µi )(λ + µm )(λ + µc ) − bp ρLs ζβLh L2
=
λ4 + aλ3 + bλ2 + cλ + d,
(23) where
a b c d
= = = =
δ + µi + µm + µc , δµi + µm µc + (δ + µi )(µm + µc ), µm µc (δ + µi ) + δµi (µm + µc ), δµi µm µc − bp ρLs ζβLh L2 .
Following a theorem due to Strelitz (1977), the polynomial f1 (λ) is stable if and only if a > 0, b > 0, 0 < c < ab, 0 < d < (abc − c2 )/a2 . Note it is clear that ¯ 0 < 1 (which implies d < 0), then EDF E is unstable. a > 0, b > 0, c > 0. If R ¯ If R0 > 1, then d > 0, we have to check the other restrictions in order to test whether f1 (λ) is stable or not. After some calculation and simplification, we have ab − c =
[(δ + µi )δµi + (δ + µi )2 (µm + µc ) + (µm + µc )µm µc +(µm + µc )2 (δ + µi )
>
0.
=
c(ab − c)
=
δµi µm µc (δ + µi )2 + (δ + µi )(µm + µc )(δµi )2
(24)
Notice that abc − c2
+µm µc (µm + µc )(δ + µi )3 + δµi (δ + µi )2 (µm + µc )2
(25)
+(µm + µc )(µm µc )2 (δ + µi ) + δµi µm µc (µm + µc )2 +µm µc (µm + µc )2 (δ + µi )2 + δµi (δ + µi )(µm + µc )3 , and a2 d = (δ + µi + µm + µc )2 [δµi µm µc − bp ρLs ζβLh L2 ] (26) +2δµi µm µc (δ + µi )(µm + µc ) − bp ρLs ζβLh L2 (δ + µi + µm + µc )2 . 14
Thus abc − c2 − a2 d = δµi µm µc (δ + µi )2 + (δ + µi )(µm + µc )(δµi )2 +µm µc (µm + µc )(δ + µi )3 + δµi (δ + µi )2 (µm + µc )2 +(µm + µc )(µm µc )2 (δ + µi ) + δµi µm µc (µm + µc )2 +µm µc (µm + µc )2 (δ + µi )2 + δµi (δ + µi )(µm + µc )3 −δµi µm µc (δ + µi )2 − δµi µm µc (µm + µc )2 −2δµi µm µc (δ + µi )(µm + µc ) +bp ρLs ζβLh L2 (δ + µi + µm + µc )2 ≥ (δ + µi )(µm + µc )(δµi )2 + δµi (δ + µi )2 (µm + µc )2 +(µm + µc )(µm µc )2 (δ + µi ) + µm µc (µm + µc )2 (δ + µi )2 +(δ + µi )(µm + µc )[µm µc (δ + µi )2 + δµi (µ2m + µ2c )] > 0.
(27)
¯ 0 < 1, the above inequalities imply that the characteristic polynomial f1 (λ) If R is stable, and so is f (λ). Then, EDF E is locally asymptotical stable. Remark 4.3. We conjecture that EDF E is actually a globally asymptotical ¯ 0 < 1. Numerical simulations verified that the reproductive steady state if R ¯ 0 > 1, number is indeed a threshold for the system. As shown in Figure 4, if R ¯ 0 < 1, then parasitism disappears. the parasitism persists, but if R
5
Model with Spatial Structure
In Section 3, we have established the existence and uniqueness of a solution of System (1). In the following we mainly focus on some stability and bifurcation analysis of the system. First of all, we focus on stability at boundary steady states in terms of the threshold (reproductive number). From System (3), at steady state, m satisfies −v
dm + bp RH − µm m = 0. dx
(28)
This can be solved exactly as m(x) =
µm bp RH (1 − e− v x ). µm
(29)
Integrating in x from 0 to L, we have M=
v − µm L bp [L + (e v − 1)]RH , ¯bp RH = ¯bp P. µm µm
Similarly, at steady state, c and C can be solved as β(T −S) − µc x c(x) = µc (1 − e v ), C=
β µc [L
+
v − µvc L µc (e
¯ − C) = βI. ¯ − 1)](T − C) , β(T 15
(30)
(31)
Mean Parasite Load / Person 0.07 0.065 0.06
average load
0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02
0
500
1000
1500
time (in year) Mean Parasite Load / Person 15
average load
10
5
0
0
5
10 time (in year)
15
20
¯ 0, R ¯0 = Figure 4: Mean parasite load as reproductive number for different R ¯ 1.000622 on the top and R0 = 0.999378 on the bottom, where Lh = 800, and r ¯ 0. are chosen to match R
16
It is evident that E0 , E1 , E2 and EDF E are boundary equilibria of System (3), too. Moreover, we have the same stability results. Theorem 5.1. E0 , E1 and E2 are unstable. Proof. By contradiction. Suppose the equilibrium E ∗ = (H ∗ , R∗ , S ∗ , T ∗ , m∗ , c∗ ) is stable, then it is a steady state of the following system
dH dt
=
(bh − µh0 − µh1 H)H − αRH,
dR dt
=
ζC −
dS dt dT dt
α 2 R − (bh + µp + α + r)R, κ
=
(bs − µs0 − µs1 T )S − ρM S,
=
(bs − µs0 − µs1 T )T − (bs + ds )(T − S),
(32)
¯ − S) and M = ¯b HR. Its Jacobian is where C = β(T p J4 (E) =
rh − 2µh1 H − αR 0 ¯ −ρbp RS 0
−αH 2α − R−δ κ −ρ¯bp HS 0
0 −ζ β¯
0 ζ β¯
rs − µs1 T − ρ¯bp HR µi
−µs1 S rs − 2µs1 T − µi
By checking the eigenvalues of J4 at E0 , E1 and E2 , we prove the conclusion. Define the reproductive number for System (1), or System (3), as R0 =
¯ ¯b L L ζ βρ p h s . δµi
(33)
Theorem 5.2. If R0 < 1, then EDF E is locally asymptotically stable; if R0 > 1 then EDF E is unstable. Proof. Refer to the previous theorem, stable equilibria satisfies the same dynamic system. The proof can be done following the same process as in [11] using the above Jacobian J4 . Concerning the interior equilibria, we have exactly same results as for the ¯ ¯bp by β¯ and ¯bp respectively. spatially independent model, replacing β, Remark 5.1. Though we have not proven it yet, EDF E seems to be globally asymptotically stable if R0 < 1.
17
.
Remark 5.2. It is worth to notice that R0 is a decreasing function with respect to velocity v and treatment rate r. This can be seen as follows, dβ¯ dv
and
=
µc d β v [ {L + (e− v L − 1)}] dv µc µc
=
v µc − µc L β 1 − µc L [ (e v − 1) + Le v ] µc µc µc v 2
=
µc β 1 L 1 [( + )e− v L − ] µc µc v µc
=
β − µc L 1 [e v − ], µc 1 + µvc L β L, lim+ β¯ = µc v→0 lim β¯ = 0.
(34)
(35)
v→∞
It is clear that for any v > 0 and L > 0, β¯ is a decreasing positive function of v. ¯ it is also a decreasing positive Noticing that ¯bp has same kind of structure as β, function of v. Then, R0 is a decreasing function of v. It is obvious that it is a decreasing function of r, too. Concerning the positive interior equilibria, the bifurcation diagram is given as follows, where parameters are chosen as in the previous section, and Lh = 800, µm = µc = 180. As we have shown, the reproductive number R0 is a decreasing function with respect to the parameters r and v. Figure 5 shows the bifurcation diagram for the system with respect to these parameters. Figure 6 is the region enclosed by dotted line in Figure 5 in a zoom-in scale. In Region I, R0 < 1, then EDF E is the only non-trivial equilibrium and numerical simulation shows that it is the only stable equilibrium. In Region II, there is one interior equilibrium which is globally stable. In region III, there are three interior equilibria, and there are two of them are asymptotically stable and no oscillatory solution exists; In region IV, there are three interior equilibria, but two adjacent ones are unstable (EDF E is also not stable); therefore, oscillatory solutions may occur in this case. In the simulations, we apply the Runge-Kutta Method for O.D.E systems and a Crank–Nicolson-type Finite Difference Method of Characteristics. It can be shown easily that the scheme is second order. Spatial dynamics of parasite larvae is determined by the river velocity v as the convection coefficient in our model, therefore, miracidia and cercariae should have same type of spatial distribution near steady states. In Figure 7 (µm = µc = 1, v = 1), we show profiles of miracidia with different initial distribution and the cercariae follow the same pattern.
18
Bifurcation Diagram I
River Current Velocity v
150
100 I
II
50
0
0
2
4
6
8 10 treatment rate r
12
14
16
Figure 5: Bifurcation diagram with respect to parameters v and r.
Bifurcation Diagram II 110 100
River Current Velocity v
90 II
80 70 60 50 40
III
30 20 10 0
IV 0
0.1
0.2
0.3
0.4 0.5 0.6 treatment rate r
0.7
0.8
0.9
1
Figure 6: Zoom-in into the dotted region in Figure 5.
19
Miracidia 5
x 10 12
density m(x,t)
10 8 6 4 2 0 1 2 1.5
0.5
1 0.5
space x
0
0
time t
Miracidia 5
x 10 12
density m(x,t)
10 8 6 4 2 0 1 2 1.5
0.5
1 0.5
space x
0
0
time t
Figure 7: Profiles of miracidia as sine (top) and cubic (bottom) initial distribution. Parameters used in the simulation are same as before, but µm = µc = 1, and v = 1.
20
6
Discussions and Conclusions
It has been observed that schistosomiasis is distributed non-uniformly in space and higher prevalence of disease is associated with low velocity of river current. Faster river currents may flush away the intermediate snail hosts thus decreasing the probability of amplification of miracidia by contacting snails. On the other hand, at the same time the river current flushes away snails, it also flushes away the parasite larvae. It is known that the miracidia and cercariae are very active free-swimming larvae. Some prefer intense light and try to stay on the surface of the river, which gives them a chance to flow with the river. So far it is not clear whether low incidence of schistosomiasis is mainly due to the flush-away of snails or of larvae. It has been observed that certain kinds of snails stay on the bed whenever they settle down, but some migrate against the river. In this paper we neglected all these movement, assuming that the distribution of snails in space is uniform. Our model clearly shows the fact that the river flow can flush out parasite larvae, making the mean parasite load decreasing with respect to river current velocity. In our model the river current velocity is assumed constant. In order to adjust to non-uniform river velocities, we have to make a change in the first four equations. One possible way is to divide the system according to different ranges of velocity and still use ordinary differential equations. Another way is to use partial differential equations which means that every species are moving along spatial direction. Numerical simulations show that it is beneficial for human hosts and snails to reside in regions of faster flow or near the (uninfected) sourcde of the river.
References [1] R.M. Anderson and R.M. May, Prevalence of schistosome infections within molluscan populations: Observed patterns and theoretical predictions, Parasitol. 79 (1979), pp. 63C94 [2] Joel E. Cohen Mathematical Models of Schistosomiasis. Ann. Rev. Ecol. Syst., 1977.8:209-233. [3] Feng, Z.; Eppert, A.; Milner, F. A.; Minchella,D.J. Estimation of parameters governing the transmission dynamics of schistosomes. Appl. Math. Lett. 17 (2004), no. 10, 1105–1112. 92D25 (93E10) [4] Z. Feng, C.-C. Li, and F. A. Milner, Schistosomiasis models with density dependence and age of infection in snail dynamics. Deterministic and stochastic modeling of biointeraction. Math. Biosci. 177/178 (2002), 271-286. [5] Feng,Z.;and Milner,F.A. A new mathematical model of schistosomiasis. Mathematical models in medical and health science (Nashville, TN, 1997), 117–128, Innov. Appl. Math., Vanderbilt Univ. Press, Nashville, TN, 1998.
21
[6] FINGERUT J.; ZIMMER C.; ZIMMER R.; Larval swimming overpowers turbulent mixing and facilitates transmission of a marine parasite. Ecology, 2003, Vol. 84, No. 9, pp. 2502C2515. [7] Schwartz E, Kozarsky P, Wilson M, Cetron M. Schistosome infection among river rafters on Omo River, Ethiopia, J Travel Med. 2005 Jan-Feb;12(1):38. [8] UTZINGER J.; MAYOMBANA C.; SMITH T.; TANNER M.; Spatial microhabitat selection by Biomphalaria pfeifferi in a small perennial river in Tanzania. Hydrobiologia (Hydrobiologia) 1997, vol. 356, pp. 53-60 [9] P. Zhang, Z. Feng, and F. A. Milner, Mathematical models of schistosomiasis with an age structure in the human host. Math. Biosc. 205 (2007), 83-107. [10] WHO TDR 2005, Seventeenth Programme Report, Progress 2003-2004. [11] Milner F.; Zhao R. A mathematical model of schistosomiasis, revisiting a model of 1997. Submitted
22
1
Introduction
Schistosomiasis, a parasite (Schistosome)-induced disease, is also known as bilharziasis after Theodor Bilharz, who first identified the parasite in Egypt in 1851. Infection is widespread with a relatively low mortality rate, but a high morbidity rate, causing severe debilitating illness. According to WHO [10], an estimated 170 million people in Sub-Saharan Africa, and a further 30 million in North Africa, Asia and South America, suffered from schistosomiasis in 2004. The disease is generally associated with rural poverty. Schistosomes have to go through an intermediate host (snails in most cases) to complete their life-cycle: eggs, miracidia, cercariae, to the final stage of adult flukes. There have been many studies in mathematical modeling using unstructured models [1, 2, 5, 11], and age-structured models [4, 9]. However, none of them incorporated spatial structure in the habitat, such as a river system. Indeed, the river current plays an important role in spreading schistosomiasis, either by affecting intermediate host snails or parasite larvae. It has been noticed in [7] that the probability for explorers in Africa to be infected in the Southern segment of the Omo River is larger than in the Northern part. Two 1
main features in the Southern Omo River are believed to contribute to the increased incidence, slow flow and large number of residences. In [8] Utzinger, et al. have observed that certain snails (such as B. pfeifferi) have preference for a certain range of river current velocities, which may explain the low incidence of schistosomiasis due to lack of snails. Fingerut and his colleagues, [6], found that environmentally triggered downward swimming can quickly bring larvae (Planktonic cercariae) to the bed, promoting contact with benthic intermediate hosts under most flow conditions, which explains the prevalence of infected snails. The main aim in this paper is to investigate how the river current can change the dynamics of schistosomiasis by considering spatial structure only in parasite larvae, namely miracidia and cercariae. The basic structure of the mathematical model is still as in [1, 5, 11], with the following assumptions: (H1) Human beings grow logistically in the absence of parasites; (H2) Snail population grows logistically in the absence of parasites; (H3) Infected Snails do not reproduce; (H4) Distribution of parasites among human hosts is overdispersed (negative binomial); (H5) Birth rates of hosts are greater than their death rates, bj −µj0 > 0, j = h, s. It is natural to incorporate spatial structure for parasite larvae first of all, since miracidia and cercariae are free-swimming in river. However, for convenience and in order to have a simpler model, human hosts, parasites and snails will be considered spatially independent. Since human beings are active more or less everywhere in their habitat, their contribution to the environment can be viewed as approximately homogeneous. The same is true for parasites since they reside inside the human body. The present paper is structured as follows: in the next section, a mathematical model with spatial structure is given; in Section 3, its well-posedness is established; in Section 4, a reduced model—spatially uniform, given by a dynamical system—is studied analytically and numerically; in Section 5, the full spatial model is partly analyzed and simulations are provided. Finally, we present some discussion and draw some conclusions in Section 6.
2
A Mathematical Model with Spatial Structure
Schistosomiasis is usually prevalent in villages near a river or a lake. For simplicity, we assume the aquatic habitat to be a river [0, L], where 0 is the origin of the river and L represents the end of the interested segment of river. Thus, miracidia and cercariae move along with the river current with velocity v. Moreover, we assume that there is no parasitism at the origin of the river, while at the end of the river, the parasitism is determined by the system. 2
In [11] it was shown that the long-term exponential growth of humans used in [5] is not adequate, not only based on demographic considerations, but also on epidemiological ones since the mean parasite load stays the same irrespective of treatment rate unless such rate is sufficient to kill all parasites. In this paper, we focus on the logistic growth in the human population, which is more adequate for long-term demographic projection and it also does away with the problem of mean parasite load not diminishing with treatment. Similar results are expected if we change the logistic model to one with constant recruitment. For simplicity, the velocity of the river flow is assumed to be constant. It has been observed that the miracidia and cercariae move in the river with directional preference—e.g. towards to light [6], and towards intermediate hosts snails. In this paper, we neglect these anisotropic movement and assume that the larvae move only due to the river movement. In other words, the spatial structure is reflected on the assumption that the parasite larvae follow a first order hyperbolic equation with convection coefficient exactly equal to the velocity of river current. Prior to building the model, let us introduce some demographic and epidemiological parameters. bh µh0 µh1 Lh α β λ κ µp r bs µs0 µs1 Ls ρ ds bp ζ µm µc v
= = = = = = = = = = = = = = = = = = = = =
per capita birth rate of human beings per capita natural death rate of human beings additional per capita death rate of human beings due to overcrowding carrying capacity of environment for human beings parasite-induced death rate of human beings per parasite production rate of cercariae per infected snail attachment rate of cercariae to human beings “clumping parameter” of binomial distribution of adult parasites per capita death rate of parasites treatment-induced rate of parasites per capita birth rate of snails per capita natural death rate of snails additional per capita death rate of snails due to overcrowding carrying capacity of environment for snails contact rate (per snail) of uninfected snails with miracidia parasite-induced death rate of snails mean number of eggs (miracidia) laid per parasite “efficacy” of cercariae in infecting humans (0 ≤ ζ ≤ 1) per capita death rate of miracidia per capita death rate of cercariae river current velocity
Let H(t), P (t), S(t), I(t) denote, respectively, the total numbers of human hosts, adult parasites, uninfected snails and infected snails at time t. Let m(x, t), c(x, t) denote, respectively, the density of free-living miracidia and cercariae at Z L Z L time t and location x, and C(t) = c(x, t)dx, M (t) = m(x, t)dx denote, 0
0
3
respectively, the total numbers of miracidia and cercariae at time t. The mathematical model we propose consists of the following system of differential equations:
dH dt
= (bh − µh0 − µh1 H)H − αP,
dP dt
=
dS dt
= (bs − µs0 − µs1 (S + I))S − ρM S,
dI dt ∂m ∂t ∂c ∂t
ζHC − α
¡ κ+1 ¢ P 2 κ
H
− (µh0 + µh1 H + µp + α + r)P,
(1) =
ρM S − (µs0 + µs1 (S + I) + ds )I,
=
−v
∂m + bp P − µm m, ∂x
=
−v
∂c + βI − µc c ∂x
along with homogeneous Dirichlet boundary conditions at the inflow boundary for the two partial differential equations, m(0, t) = 0,
c(0, t) = 0.
(2)
Remark 2.1. The singularity at H = 0 in System (1) can be avoided by introducing two new variables, the mean parasite load per person R = P/H, and total population size of snails T = S + I. It is easy to derive the following equivalent system to System (1) describing the evolution of H, R, S, T, m, c:
dH dt
=
(bh − µh0 − µh1 H)H − αRH,
dR dt
=
ζC −
dS dt
=
(bs − µs0 − µs1 T )S − ρM S,
dT dt ∂m ∂t ∂c ∂t
α 2 R − (bh + µp + α + r)R, κ
(3) =
(bs − µs0 − µs1 T )T − (bs + ds )(T − S)
=
−v
∂m + bp RH − µm m, ∂x
=
−v
∂c + β(T − S) − µc c. ∂x 4
Remark 2.2. Free-living miracidia and cercariae can live up to 48 hours in fresh water, and they are actively searching hosts only in the first several hours. Therefore, we can neglect the small portion of those penetrating into hosts and assume the loss of larvae is only due to natural death. Remark 2.3. Homogeneous Dirichlet boundary conditions in this model are assumed for simplicity and they can be replaced by measured field data when a specific region needs to be considered.
3
Well-Posedness
Concerning the well-posedness of System (1), we can equivalently investigate System (3). Since (3.v) and (3.vi) are first order hyperbolic equations, we can solve them exactly along characteristic lines x = τ + vt, where τ is a parameter. Denote m(t) ˜ = m(τ + vt, t); then Eq. (3.v) is equivalent to dm ˜ + µm m ˜ = bp RH, dt
(4)
and also to
d µm t (e m) ˜ = bp eµm t RH. (5) dt Since the habitat has finite length L, we can group characteristic lines as above x = vt and below x = vt. In the region {(x, t) : vt < x < L}, integrating Eq. (5) from 0 to t, we have Z µm t
e
t
m(t) ˜ = m(0) ˜ + bp
eµm s R(s)H(s)ds,
(6)
0
and then Z
t
m(x, t) = e−µm t m(x − vt, 0) + bp
e−µm (t−s) R(s)H(s)ds.
(7)
0
x Similarly, in the region {(x, t) : 0 < x < vt}, integrating Eq. (5) from t − to v t, we have Z t x x eµm t m(t) ˜ = eµm (t− v ) m(t ˜ − ) + bp eµm s R(s)H(s)ds, (8) x v t− v and then, m(x, t) =
− µvm x
e
Z =
bp 0
x v
x m(0, t − ) + bp v −µm ( x v −s)
e
Z
t
e−µm (t−s) R(s)H(s)ds
t− x v
x x R(s + t − )H(s + t − )ds. v v 5
(9)
In summary, we can write the solution as Z t −µm t e m(x − vt, 0) + b e−µm (t−s) R(s)H(s)ds if x > vt, p 0 m(x, t) = Z xv x x x bp if x ≤ vt. e−µm ( v −s) R(s + t − )H(s + t − )ds v v 0
(10)
Integrating Eq. (10) in the spatial variable x on the interval [0, L], the total Z L population of miracidia M (t) = m(x, t)dx at time t is given by 0
Z vt Z xv x x x e−µm ( v −s) R(s + t − )H(s + t − )dsdx if t < b p v v 0 0 ¸ Z t Z L· + e−µm t m(x − vt, 0) + bp e−µm (t−s) R(s)H(s)ds dx, vt 0 x Z Z L v x x x bp e−µm ( v −s) R(s + t − )H(s + t − )dsdx if t ≥ v v 0 0
L , v (11) L . v
Similarly, we can write the solution of Eq. (3.vi) as Z t −µc t c(x − vt, 0) + β e−µc (t−s) (T (s) − S(s))ds if x > vt, e 0 Z xv c(x, t) = x x x β e−µc ( v −s) (T (s + t − ) − S(s + t − ))ds if x ≤ vt, v v 0 Z L and the total population of cercariae C(t) = c(x, t)dx at time t as
(12)
0
Z L· ¸ Z t −µc t −µc (t−s) e m(x − vt, 0) + β e (T (s) − S(s))ds dx vt 0 Z vt Z xv x x x e−µc ( v −s) (T (s + t − ) − S(s + t − ))dsdx if t < +β v v 0 0 Z L Z xv x x x β e−µc ( v −s) (T (s + t − ) − S(s + t − ))dsdx if t ≥ v v 0 0
L v,
(13)
L v.
Concerning the positivity of solutions to System (3), subtracting Eq. (3.iii) from Eq. (3.iv), we have d(T − S) = (bs − µs0 − µs1 T )(T − S) − (bs + ds )(T − S) + ρM. dt
(14)
System (3) implies that H and S are nonnegative if the initial data are nonnegative. It is clear from Eq. (14) that if M is nonnegative, then T − S 6
is nonnegative. Moreover, from Eq. (3), if T − S is nonnegative, then R is nonnegative. Eq. (11) shows that M is nonnegative if R is nonnegative, and Eq. (13) shows that C is nonnegative if T − S is nonnegative. Therefore, when H(0) > 0 and S(0) > 0, if either R(0) > 0 or T (0) > S(0), then for small ε > 0 we certainly have H, R, S, M, C and T −S positive on (0, ε). It is clear that Eq. (3.i) gives the positivity of H, Eq. (3.iii) gives the positivity of S, and then Eq. (3.iv) gives the positivity of T as long as they exist. Hence, it suffices to show that T − S must remain positive as long as it exists, since that is sufficient (in fact, equivalent) to establishing the positivity of R. Suppose that T − S vanishes for some positive time, and let t¯ be the smallest such time. Then, T − S is positive on [0, t¯) and T (t¯) = S(t¯). It follows that C(t¯) > 0 from Eq. (13). Then, R is positive on [0, t¯] from Eq. (3.ii) and then Eq. (14) implies that T − S is positive on [0, t¯], contradicting T (t¯) = S(t¯). Thus T − S remains positive as long as it exists, and so do R, M and C. Moreover, m(x, t) and c(x, t) are also positive if R(t), H(t), T (t) − S(t) are positive. This completes the proof of positivity of solutions to System (3). Concerning their boundedness, it is evident that Eq. (3.a) and Eq. (3.d) implies that H ≤ Lh and T ≤ Ls . Moreover, S ≤ Ls and I ≤ Ls . For convenience, let us introduce a notation δ = bh + µp + α + r. It is not hard to show that C is bounded by observing that for t > L/v, Z |C(t)| =
L
Z
x v
β 0
x
e−µc ( v −s) (T (s + t −
0
Z
L
≤
βLs µc
=
· ¸ µc βLs v L + (e− v L − 1) µc µc
,
¯ C.
x x ) − S(s + t − ))dsdx v v
x
(1 − e−µc v )dx (15)
0
Notice from Eq. (3.ii) that α dR = ζC − R2 − δR ≤ ζ C¯ − δR, dt κ
(16)
The boundedness of R follows from R ≤ R0 e−δt +
ζ C¯ ζ C¯ ¯ (1 − e−δt ) ≤ R0 + , R. δ δ
7
(17)
Last, for the boundedness of M , it is sufficient to check that for t > L Ã /v, Z L Z xv x x x |M (t)| = bp e−µm ( v −s) R(s + t − )H(s + t − ))dsdx v v 0 0 ≤
¯Z L x bp Lh R (1 − e−µm v )dx µm 0
=
¸ ¯· bp Lh R v − µc L L+ (e v − 1) µm µm
,
¯. M
(18)
Now let us introduce a vector V¯ = (H, R, S, T )t . Then System (3) can be rewritten as dV¯ = F¯ (V¯ ) = (f1 , f2 , f3 , f4 )t , (19) dt where f1 (H, R, S, T ) =
(bh − µh0 − µh1 H)H − αRH,
f2 (H, R, S, T ) =
ζC −
α 2 R − (bh + µp + α + r)R, κ
f3 (H, R, S, T ) = (bs − µs0 − µs1 T )S − ρM S, f4 (H, R, S, T ) =
(bs − µs0 − µs1 T )T − (bs + ds )(T − S),
and M and C are functions of T − S, H, and R defined in (11) and (13). It is easy to show that F¯ is Lipschitz continuous and Picard Iteration Theorem (fixed-point theory) proves local existence and uniqueness of the solution. The global existence and uniqueness of solutions is guaranteed by their a priori boundedness. Then, the following theorem has been proved. Theorem 3.1. System (1), or equivalently, System (3) admits a unique solution for all time. Moreover, if H(0) and S(0) are positive and either T (0) > S(0), or M (0) > 0, or R(0) > 0, or C(0) > 0, then H, R, S, m and c are positive and T > S as long as they exist.
4
Special Case: Model in Spatially Uniform Distribution
At the first glance, it seems strange to talk about a case where all species are uniformly distributed in space since the homogeneous in-flow boundary conditions for miracidia and cercariae will force a disease free solution. However, spatially uniform distribution can happen at steady state if non-homogeneous in-flow boundary conditions are imposed. 8
The diversity of steady states for the uniformly distributed model has been established in [11], except for β in [11] replaced by βL/µc and bp by bp L/µm . In[11] the population size of miracidia is assumed to be simply proportional to that of parasites and the population size of cercariae is proportional to that of infected snails. It is clear that System (3) also admits four parasite-free steady states E0 E1 E2 EDF E
= (H0 , R0 , S0 , T0 , m0 , c0 ) = (H1 , R1 , S1 , T1 , m1 , c1 ) = (H2 , R2 , S2 , T2 , m2 , c2 ) = (HDF E , RDF E , SDF E , TDF E , mDF E , cDF E )
= (0, 0, 0, 0, 0, 0), = (0, 0, Ls , Ls , 0, 0), = (Lh , 0, 0, 0, 0, 0), = (Lh , 0, Ls , Ls , 0, 0).
where, EDF E represents (H, P, S, I, m, c) = (Lh , 0, Ls , 0, 0, 0) for System (1). E0 is the trivial equilibrium, while E1 and E2 are semi-trivial in that they represent, respectively, the ultimate state of a logistic population of snails and of humans in the absence of parasitism. Concerning the non-trivial equilibria, as in [11], we see from System (3) that at steady state H satisfies the cubic equation H 3 + a2 H 2 + a1 H + a0 = 0, where
a2
=
−Lh ,
a1
=
rs α αµs1 µi + ¯ ¯bp )2 , ρ¯bp µh1 κζ β(ρ
a0
=
−
(20)
δαµs1 µi rh αµs1 µi ¯ h1 (ρ¯bp )2 − ζ βµ ¯ h1 (ρ¯bp )2 . κζ βµ
¯ µi are defined as Here, ¯bp , β, ¯bp = bp L , β¯ = βL , µi = bs + ds . µm µc Define the discriminant of the above cubic polynomial as D = Q3 + R2 , where Q=
3a1 − a22 , 9
U=
9a1 a2 − 27a0 − 2a32 . 54
Note that a2 < 0, a1 > 0, and a0 < 0, and we have a similar result as in [11]. Theorem 4.1. Let D, Q and U be defined as above, and consider D and U as functions of treatment rate r. • If Q > 0 or if Q < 0, U (0) > 0 and D(0) > 0, then Eq. (20) admits only one positive root; • If Q < 0, and if D(0) < 0, then Eq. (20) admits three, then two, and finally one positive root as r increases. 9
• If Q < 0, U (0) < 0 and D(0) > 0, then Eq. (20) admits one, two, then three, two, and finally one positive root as r increases. Remark 4.1. It can be easily seen fromµSystem ¶ (3) that if a root H of Eq. (20) r H ρ¯bp rh Lh h satisfies 0 < H < Lh and rs > ρ¯bp 1− H (or rs > ), then α Lh 4α System (3) or, equivalently, System (1) has a positive equilibrium with human population size H.
Bifurcation Diagram 1200
Carrying capacity of human Lh
1000 IV 800 III 600
400
II
200
0
I 0
0.2
0.4 0.6 treatment rate r
0.8
1
Figure 1: Bifurcation Diagram. In region I, there is no interior positive equilibria and EDF E is globally asymptotically stable. In region II, there is one and only one interior positive equilibrium and it is globally asymptotically stable. In region III, there are three positive interior equilibria and two of them are locally asymptotically stable and no oscillatory solution happens. In region IV, there are also three positive interior equilibria and two of them are unstable and oscillatory solution occurs for certain initial data. Figure 1 shows the diversity of system (1) in the number of positive interior equilibria with respect to two bifurcation parameters Lh and r. The parameters are chosen as: L = 1, µh0 = 1/70, bh = 3/140, α = 0.00001, ζ = 0.000027, κ = 0.243, µp = 0.2, µs0 = 0.5, bs = 80, Ls = 20000, ρ = 0.00004, ds = 6.0, µm = 180, µc = 180, β = 4000µc , bp = 20µm . Remark 4.2. Indeed, Figure 1 is a bifurcation diagram of System (1) for the 10
5
Human Hosts
population size
1.2
625
624.5
624
0 6
population size
8
x 10
500 time (in year) Larvae
4 2
0
Parasites
1 0.9
0
500 time (in year) Snails
1000
2000 M C
6
0
x 10
1.1
0.8
1000
population size
population size
625.5
500 time (in year)
1000 500 0
1000
S I
1500
0
500 time (in year)
1000
Figure 2: Oscillatory solution when r = 0 and Lh = 800 on [0, 1000].
5
Human Hosts
population size
1.2
625.45
625.4
625.35 990 6
population size
8
x 10
995 time (in year) Larvae
Parasites
1.1 1 0.9
995 time (in year) Snails
1000
2000 M C
6 4 2 0 990
x 10
0.8 990
1000
population size
population size
625.5
995 time (in year)
1000
S I
1500 1000 500 0 990
995 time (in year)
1000
Figure 3: Oscillatory solution when r = 0 and Lh = 800 on [990, 1000].
11
chosen parameters. We have been unable to derive a stability analysis for positive equilibria. Numerical simulations show that System (1) may converge to a steady state or oscillatory solutions. Figure 2 and Figure 3 show that an oscillatory solution occurs when parameters are in region IV, where the simulations were done using r = 0 and Lh = 800. However, concerning the stability of boundary steady states, we obtain here similar results as in [11], too. Theorem 4.2. E0 , E1 and E2 are unstable. Proof. For convenience, introduce another notation: κ+1 . κ The Jacobian of System (3) at the equilibrium E = (H, R, S, T, m, c) is a11 −αH 0 0 0 0 0 a22 0 0 0 ζL 0 0 a −µ S −ρLS 0 33 s1 . J2 (E) = 0 (bs + ds ) a44 0 0 0 b p R bp H 0 0 −µm 0 0 0 −β β 0 −µc α ¯=α
where
Then,
J2 (E0 ) =
rh 0 0 0 0 0
a11
=
rh − 2µh1 H − αR,
a22
=
−
a33
=
rs − µs1 T − ρLm,
a44
=
rs − 2µs1 T − (bs + ds ).
0 −δ 0 0 0 0
2α R − δ, κ
0 0 rs (bs + ds ) 0 −β
0 0 0 0 0 0 −(µs0 + ds ) 0 0 −µm β 0
0 ζL 0 0 0 −µc
.
It is obvious that λ = rh > 0 is an eigenvalue of J3 (E0 ); then E0 is unstable. Next, rh 0 0 0 0 0 0 −δ 0 0 0 ζL 0 0 0 −µ L −ρLL 0 s1 s s . J2 (E1 ) = 0 (bs + ds ) −rs − (bs + ds ) 0 0 0 0 0 0 0 −µm 0 0 0 −β β 0 −µc 12
Again λ = rh > 0 is an eigenvalue of J3 (E1 ) and then E1 is unstable. Finally, −rh −αLh 0 0 0 0 0 −δ 0 0 0 ζL 0 0 rs 0 0 0 J2 (E2 ) = 0 0 (b + d ) −(µ + d ) 0 0 s s s0 s 0 b p Lh 0 0 −µm 0 0 0 −β β 0 −µc
.
Note that λ = rs > 0 is an eigenvalue of J3 (E2 ) and then E2 is unstable. Regarding the stability of the Disease Free Equilibrium, we have the following result. Define so-called reproductive number of parasitism ¯ ¯ ¯ 0 = ζ βρbp Lh Ls . R (21) δµi ¯ 0 < 1, then EDF E is locally asymptotically stable; if R ¯0 > 1 Theorem 4.3. If R then EDF E is unstable. Proof. Note that System (1) and System (3) are equivalent. It turns out to be easier to find the threshold in this case by looking at the Jacobian of System (1) than System (3). The Jacobian of System (1) at the equilibrium (H, P, S, I, m, c) is rh − 2µh1 H −α 0 0 0 0 a21 a22 0 0 0 ζLH 0 0 a −µ S −ρLS 0 33 s1 , J1 (E) = 0 0 ρLm − µs1 I a44 ρLS 0 0 bp 0 0 −µm 0 0 0 0 β 0 −µc where
µ a21
= ζLc + α ¯ µ
2P H
P H
¶2 − µh1 P,
¶
a22
= −¯ α
− (µh0 + µh1 H + µp + α + r),
a33
= rs − µs1 (2S + I) − ρLm,
a44 = −(µs0 + µs1 (S + 2I) + ds ). At the equilibrium EDEF , the Jacobian is −rh −α 0 0 0 −δ 0 0 0 0 −r −r s s J1 (EDF E ) = 0 0 0 −(b + ds ) s 0 bp 0 0 0 0 0 β 13
0 0 −ρLLs ρLLs −µm 0
0 ζLLh 0 0 0 −µc
.
The characteristic polynomial of the above Jacobian is f (λ) = (λ + (bh − µh ))(λ + (bs − µs )) (22) [(λ + δ)(λ + µi )(λ + µm )(λ + µc ) − bp ρLs ζβLh L2 ]. It is clear that λ = −(bh − µh ) and λ = −(bs − µs ) are two roots of f (λ) and they are negative from our assumptions. Whether f (λ) has positive roots is completely determined by the forth order multiplier f1 (λ)
=
(λ + δ)(λ + µi )(λ + µm )(λ + µc ) − bp ρLs ζβLh L2
=
λ4 + aλ3 + bλ2 + cλ + d,
(23) where
a b c d
= = = =
δ + µi + µm + µc , δµi + µm µc + (δ + µi )(µm + µc ), µm µc (δ + µi ) + δµi (µm + µc ), δµi µm µc − bp ρLs ζβLh L2 .
Following a theorem due to Strelitz (1977), the polynomial f1 (λ) is stable if and only if a > 0, b > 0, 0 < c < ab, 0 < d < (abc − c2 )/a2 . Note it is clear that ¯ 0 < 1 (which implies d < 0), then EDF E is unstable. a > 0, b > 0, c > 0. If R ¯ If R0 > 1, then d > 0, we have to check the other restrictions in order to test whether f1 (λ) is stable or not. After some calculation and simplification, we have ab − c =
[(δ + µi )δµi + (δ + µi )2 (µm + µc ) + (µm + µc )µm µc +(µm + µc )2 (δ + µi )
>
0.
=
c(ab − c)
=
δµi µm µc (δ + µi )2 + (δ + µi )(µm + µc )(δµi )2
(24)
Notice that abc − c2
+µm µc (µm + µc )(δ + µi )3 + δµi (δ + µi )2 (µm + µc )2
(25)
+(µm + µc )(µm µc )2 (δ + µi ) + δµi µm µc (µm + µc )2 +µm µc (µm + µc )2 (δ + µi )2 + δµi (δ + µi )(µm + µc )3 , and a2 d = (δ + µi + µm + µc )2 [δµi µm µc − bp ρLs ζβLh L2 ] (26) +2δµi µm µc (δ + µi )(µm + µc ) − bp ρLs ζβLh L2 (δ + µi + µm + µc )2 . 14
Thus abc − c2 − a2 d = δµi µm µc (δ + µi )2 + (δ + µi )(µm + µc )(δµi )2 +µm µc (µm + µc )(δ + µi )3 + δµi (δ + µi )2 (µm + µc )2 +(µm + µc )(µm µc )2 (δ + µi ) + δµi µm µc (µm + µc )2 +µm µc (µm + µc )2 (δ + µi )2 + δµi (δ + µi )(µm + µc )3 −δµi µm µc (δ + µi )2 − δµi µm µc (µm + µc )2 −2δµi µm µc (δ + µi )(µm + µc ) +bp ρLs ζβLh L2 (δ + µi + µm + µc )2 ≥ (δ + µi )(µm + µc )(δµi )2 + δµi (δ + µi )2 (µm + µc )2 +(µm + µc )(µm µc )2 (δ + µi ) + µm µc (µm + µc )2 (δ + µi )2 +(δ + µi )(µm + µc )[µm µc (δ + µi )2 + δµi (µ2m + µ2c )] > 0.
(27)
¯ 0 < 1, the above inequalities imply that the characteristic polynomial f1 (λ) If R is stable, and so is f (λ). Then, EDF E is locally asymptotical stable. Remark 4.3. We conjecture that EDF E is actually a globally asymptotical ¯ 0 < 1. Numerical simulations verified that the reproductive steady state if R ¯ 0 > 1, number is indeed a threshold for the system. As shown in Figure 4, if R ¯ 0 < 1, then parasitism disappears. the parasitism persists, but if R
5
Model with Spatial Structure
In Section 3, we have established the existence and uniqueness of a solution of System (1). In the following we mainly focus on some stability and bifurcation analysis of the system. First of all, we focus on stability at boundary steady states in terms of the threshold (reproductive number). From System (3), at steady state, m satisfies −v
dm + bp RH − µm m = 0. dx
(28)
This can be solved exactly as m(x) =
µm bp RH (1 − e− v x ). µm
(29)
Integrating in x from 0 to L, we have M=
v − µm L bp [L + (e v − 1)]RH , ¯bp RH = ¯bp P. µm µm
Similarly, at steady state, c and C can be solved as β(T −S) − µc x c(x) = µc (1 − e v ), C=
β µc [L
+
v − µvc L µc (e
¯ − C) = βI. ¯ − 1)](T − C) , β(T 15
(30)
(31)
Mean Parasite Load / Person 0.07 0.065 0.06
average load
0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02
0
500
1000
1500
time (in year) Mean Parasite Load / Person 15
average load
10
5
0
0
5
10 time (in year)
15
20
¯ 0, R ¯0 = Figure 4: Mean parasite load as reproductive number for different R ¯ 1.000622 on the top and R0 = 0.999378 on the bottom, where Lh = 800, and r ¯ 0. are chosen to match R
16
It is evident that E0 , E1 , E2 and EDF E are boundary equilibria of System (3), too. Moreover, we have the same stability results. Theorem 5.1. E0 , E1 and E2 are unstable. Proof. By contradiction. Suppose the equilibrium E ∗ = (H ∗ , R∗ , S ∗ , T ∗ , m∗ , c∗ ) is stable, then it is a steady state of the following system
dH dt
=
(bh − µh0 − µh1 H)H − αRH,
dR dt
=
ζC −
dS dt dT dt
α 2 R − (bh + µp + α + r)R, κ
=
(bs − µs0 − µs1 T )S − ρM S,
=
(bs − µs0 − µs1 T )T − (bs + ds )(T − S),
(32)
¯ − S) and M = ¯b HR. Its Jacobian is where C = β(T p J4 (E) =
rh − 2µh1 H − αR 0 ¯ −ρbp RS 0
−αH 2α − R−δ κ −ρ¯bp HS 0
0 −ζ β¯
0 ζ β¯
rs − µs1 T − ρ¯bp HR µi
−µs1 S rs − 2µs1 T − µi
By checking the eigenvalues of J4 at E0 , E1 and E2 , we prove the conclusion. Define the reproductive number for System (1), or System (3), as R0 =
¯ ¯b L L ζ βρ p h s . δµi
(33)
Theorem 5.2. If R0 < 1, then EDF E is locally asymptotically stable; if R0 > 1 then EDF E is unstable. Proof. Refer to the previous theorem, stable equilibria satisfies the same dynamic system. The proof can be done following the same process as in [11] using the above Jacobian J4 . Concerning the interior equilibria, we have exactly same results as for the ¯ ¯bp by β¯ and ¯bp respectively. spatially independent model, replacing β, Remark 5.1. Though we have not proven it yet, EDF E seems to be globally asymptotically stable if R0 < 1.
17
.
Remark 5.2. It is worth to notice that R0 is a decreasing function with respect to velocity v and treatment rate r. This can be seen as follows, dβ¯ dv
and
=
µc d β v [ {L + (e− v L − 1)}] dv µc µc
=
v µc − µc L β 1 − µc L [ (e v − 1) + Le v ] µc µc µc v 2
=
µc β 1 L 1 [( + )e− v L − ] µc µc v µc
=
β − µc L 1 [e v − ], µc 1 + µvc L β L, lim+ β¯ = µc v→0 lim β¯ = 0.
(34)
(35)
v→∞
It is clear that for any v > 0 and L > 0, β¯ is a decreasing positive function of v. ¯ it is also a decreasing positive Noticing that ¯bp has same kind of structure as β, function of v. Then, R0 is a decreasing function of v. It is obvious that it is a decreasing function of r, too. Concerning the positive interior equilibria, the bifurcation diagram is given as follows, where parameters are chosen as in the previous section, and Lh = 800, µm = µc = 180. As we have shown, the reproductive number R0 is a decreasing function with respect to the parameters r and v. Figure 5 shows the bifurcation diagram for the system with respect to these parameters. Figure 6 is the region enclosed by dotted line in Figure 5 in a zoom-in scale. In Region I, R0 < 1, then EDF E is the only non-trivial equilibrium and numerical simulation shows that it is the only stable equilibrium. In Region II, there is one interior equilibrium which is globally stable. In region III, there are three interior equilibria, and there are two of them are asymptotically stable and no oscillatory solution exists; In region IV, there are three interior equilibria, but two adjacent ones are unstable (EDF E is also not stable); therefore, oscillatory solutions may occur in this case. In the simulations, we apply the Runge-Kutta Method for O.D.E systems and a Crank–Nicolson-type Finite Difference Method of Characteristics. It can be shown easily that the scheme is second order. Spatial dynamics of parasite larvae is determined by the river velocity v as the convection coefficient in our model, therefore, miracidia and cercariae should have same type of spatial distribution near steady states. In Figure 7 (µm = µc = 1, v = 1), we show profiles of miracidia with different initial distribution and the cercariae follow the same pattern.
18
Bifurcation Diagram I
River Current Velocity v
150
100 I
II
50
0
0
2
4
6
8 10 treatment rate r
12
14
16
Figure 5: Bifurcation diagram with respect to parameters v and r.
Bifurcation Diagram II 110 100
River Current Velocity v
90 II
80 70 60 50 40
III
30 20 10 0
IV 0
0.1
0.2
0.3
0.4 0.5 0.6 treatment rate r
0.7
0.8
0.9
1
Figure 6: Zoom-in into the dotted region in Figure 5.
19
Miracidia 5
x 10 12
density m(x,t)
10 8 6 4 2 0 1 2 1.5
0.5
1 0.5
space x
0
0
time t
Miracidia 5
x 10 12
density m(x,t)
10 8 6 4 2 0 1 2 1.5
0.5
1 0.5
space x
0
0
time t
Figure 7: Profiles of miracidia as sine (top) and cubic (bottom) initial distribution. Parameters used in the simulation are same as before, but µm = µc = 1, and v = 1.
20
6
Discussions and Conclusions
It has been observed that schistosomiasis is distributed non-uniformly in space and higher prevalence of disease is associated with low velocity of river current. Faster river currents may flush away the intermediate snail hosts thus decreasing the probability of amplification of miracidia by contacting snails. On the other hand, at the same time the river current flushes away snails, it also flushes away the parasite larvae. It is known that the miracidia and cercariae are very active free-swimming larvae. Some prefer intense light and try to stay on the surface of the river, which gives them a chance to flow with the river. So far it is not clear whether low incidence of schistosomiasis is mainly due to the flush-away of snails or of larvae. It has been observed that certain kinds of snails stay on the bed whenever they settle down, but some migrate against the river. In this paper we neglected all these movement, assuming that the distribution of snails in space is uniform. Our model clearly shows the fact that the river flow can flush out parasite larvae, making the mean parasite load decreasing with respect to river current velocity. In our model the river current velocity is assumed constant. In order to adjust to non-uniform river velocities, we have to make a change in the first four equations. One possible way is to divide the system according to different ranges of velocity and still use ordinary differential equations. Another way is to use partial differential equations which means that every species are moving along spatial direction. Numerical simulations show that it is beneficial for human hosts and snails to reside in regions of faster flow or near the (uninfected) sourcde of the river.
References [1] R.M. Anderson and R.M. May, Prevalence of schistosome infections within molluscan populations: Observed patterns and theoretical predictions, Parasitol. 79 (1979), pp. 63C94 [2] Joel E. Cohen Mathematical Models of Schistosomiasis. Ann. Rev. Ecol. Syst., 1977.8:209-233. [3] Feng, Z.; Eppert, A.; Milner, F. A.; Minchella,D.J. Estimation of parameters governing the transmission dynamics of schistosomes. Appl. Math. Lett. 17 (2004), no. 10, 1105–1112. 92D25 (93E10) [4] Z. Feng, C.-C. Li, and F. A. Milner, Schistosomiasis models with density dependence and age of infection in snail dynamics. Deterministic and stochastic modeling of biointeraction. Math. Biosci. 177/178 (2002), 271-286. [5] Feng,Z.;and Milner,F.A. A new mathematical model of schistosomiasis. Mathematical models in medical and health science (Nashville, TN, 1997), 117–128, Innov. Appl. Math., Vanderbilt Univ. Press, Nashville, TN, 1998.
21
[6] FINGERUT J.; ZIMMER C.; ZIMMER R.; Larval swimming overpowers turbulent mixing and facilitates transmission of a marine parasite. Ecology, 2003, Vol. 84, No. 9, pp. 2502C2515. [7] Schwartz E, Kozarsky P, Wilson M, Cetron M. Schistosome infection among river rafters on Omo River, Ethiopia, J Travel Med. 2005 Jan-Feb;12(1):38. [8] UTZINGER J.; MAYOMBANA C.; SMITH T.; TANNER M.; Spatial microhabitat selection by Biomphalaria pfeifferi in a small perennial river in Tanzania. Hydrobiologia (Hydrobiologia) 1997, vol. 356, pp. 53-60 [9] P. Zhang, Z. Feng, and F. A. Milner, Mathematical models of schistosomiasis with an age structure in the human host. Math. Biosc. 205 (2007), 83-107. [10] WHO TDR 2005, Seventeenth Programme Report, Progress 2003-2004. [11] Milner F.; Zhao R. A mathematical model of schistosomiasis, revisiting a model of 1997. Submitted
22
