Schistosomiasis models with two migrating human groups

Available online at www.sciencedirect.com SCIENCE

ELSEVIER

~____~D, REC:T °

MATHEMATICAL AND COMPUTER MODELLING

Mathematical and Computer Modelling 41 (2005) 1213-1230 www.elsevier.com/locate/mcm

Schistosomiasis M o d e l s with Two Migrating H u m a n Groups Z. FENG Department of Mathematics, Purdue University 150 N. University Street, West Lafayette, IN 47907-2067, U.S.A. z f eng©math, purdue, edu

C . - C . LI Holistic Education Center, St. John's and St. Mary's Institute of Technology 499, Sec. 4, Tam King Road Tamsui, Taipei, Taiwan, R.O.C. eli©mail, sj smit. edu. tw

F . A. MILNER Department of Mathematics, Purdue University 150 N. University Street, West Lafayette, IN 47907-2067, U.S.A.

milner©purdue, edu

(Received October 2003; revised and accepted October 2003) A b s t r a c t - - W e propose, in this paper, models of schistosomiasis that incorporate several realistic features including two human habitats, migration between these, negative binomial distribution of sehistosomes within human hosts, disease-induced mortality in both human and snail hosts, and others. The qualitative and quantitative mathematical properties are analyzed. Numerical simulations help examine the dynamics and suggest some properties of these models that we were unable to prove mathematically. Partial results are extended to models with multiple human groups, and numerical simulations also support the results. Explicit thresholds for the survival of schistosomes are established. Control strategies derived from these thresholds are also discussed. @ 2005 Elsevier Ltd. All rights reserved.

1. I N T R O D U C T I O N In the developing countries, schistosomiasis is frequently a serious health problem. There are several studies describing the dynamics of schistosomiasis and other parasitic diseases in human populations, e.g., [1-8]. In [9], we described some models that incorporate human hosts, adult parasites, uninfected and infected snails, free-living miracidia, and free-living cercariae. So far, the models have focused on the interactions between one group of human hosts and schistosomes in a contaminated water resource. However, in realistic situations, there might be several human groups sharing the contaminated water resource. Take the Yang-Tzi River as an example. Being the longest river in China, there is a serious schistosomiasis problem and there are also many villages, towns, and cities along the contaminated region. This research was supported in part by NSF Grant DMS-9974389. 0895-7177/05/$ - see front matter @ 2005 Elsevier Ltd. All rights reserved. doi: 10.1016/j.mcm.2004.10.023

Typeset by AA.IS-TEX

1214

Z. FENG et al.

t m21 Village 1

Village 2 m12

Cercariae

Cercariae

Miracidia

Miracidia

Snails

WaterResource Figure 1. While modeling any ecosystem in some detail, one should begin with a simple model. We consider, in this paper, two human habitats sharing the same contaminated water resource, as schematically represented in Figure 1. In this paper, we shall show some realistic models describing the disease dynamics involving two migrating human groups. We shall also study the mathematical properties of these systems. We consider a body of water where the intermediate snail hosts live that is either stagnant or moving slowly; thus, we need not consider the movement of the snails. We shall establish an explicit threshold below which the parasites die out, and above which the disease persists. We shall establish models with multiple human groups and find some structurally similarities between the models involving two human groups and those involving n groups, though some proofs are incomplete. We shall use numerical simulations to support the mathematical results and to support some suggestions of properties that we have been unable to prove through mathematical analysis. We finally discuss some control strategies by observing the roles of the different parameters in this threshold. 2. T H E

MODELS

WITH

TWO

HUMAN

GROUPS

Let N1, N2, P1, P2, S, I, C, M denote the numbers of human hosts living in Village 1, human hosts living in Village 2, adult parasites that are hosted by human hosts in Village 1, adult parasites that are hosted by human hosts in Village 2, uninfected snails, infected snails, free-living cereariae and miracidia, respectively. Let x(t, .) denote the infection-age density of infected snails at time t, i.e., fb x(t,T)dT denotes the number of snails that at time t have been infected with miracidia for a length of time ~- E (a, b), satisfying l i m ~ _ ~ x(t, ~-) = 0. The following parameters will be used in our models (i = 1, 2), all of them positive, except c~, ds >_ 0: Ah,i is bp is As is ph is #p is #s is a is ds is p is ~ is r is

the the the the the the the the the the the

recruitment rate of human hosts in Village i; per capita birth (egg-laying) rate of adult parasites; recruitment rate of snails; per capita natural death rate of human hosts; per capita death rate of adult parasites; per capita natural death rate of snails; disease-induced death rate of human hosts per parasite; disease-induced death rate of snails; per capita (successful) rate of infection of snails by one miracidium; per capita rate of infection of human by one cercaria in Village i; releasing rate of cercariae by one snail;

Schistosomiasis Models

1215

cri is the treatment rate of human hosts in Village i; rni,j is the immigration rate of human hosts from Village i to Village j. In [9], we had the following system

d - - N = Ah -- # h N - aP, dt d p= zCN_ dt

(#h + pp + a + a ) p _ a d ~ S = b(S,I)

-

( k-k- l )

(P--~-N)

# ~ S - pMS,

(1)

~---~x(t, T) ÷ ~---~x(t, T) -- --(#~ ÷ d~)x(t, T), x(t, O) = pMS, C(t) =

/o

r(T)x(t, T) dT,

x(O, m) = Xo (T), I(t) =

/o

x(t, T) dT,

M = b~P.

We consider now the situation of two neighboring villages sharing a water resource which is contaminated by schistosomes. We assume that there is migration between these two villages. Since we want to examine long-term dynamics, it would be unrealistic to use models of exponential growth for the human and snail hosts. Therefore, we use a simple recruitment model that asymptotically stabilizes those populations at a number equal to the ratio of recruitment rate to mortality rate. Then, the dynamics for the number of human hosts in Village 1, N1, and human hosts in Village 2, N2, are given by the following system:

d

~ N 1 = Ah,1 -- phN1 -- aP1 ÷ m2,1N2

-- m l , 2 N 1 ,

(2) d N 2 = Ah,2 - #hN2 -- aP2 + ml,2N1

-- TD,2,1N 2 .

The dynamics of the adult parasites in the two villages are given by a system similar to that in (1). We shall assume that the parasites are overdispersed and they have negative binomial distributions among human hosts with clumping parameters ki. We also assume the releasing rate of cercariae is infection-age independent, i.e., r ( r ) -- r. Thus, C(t) = rI(t). Therefore, the dynamics for the adult parasites hosted by the humans living in Village 1 and Village 2, respectively, P1 and P2, are given by

-

)'

<:<

(3)

P =9 "SN2-(Vh

)

Note that, based on [9], the number of miraeidia, M, should be bp(P1 + P2) in this model. We integrate the partial differential equation for x in (1) with respect to ~- from 0 to oc, and thus, obtain an ordinary differential equation for the total number of infected snails, I d

-2-I = pbp(P1 + P2)S - (Ms ÷ ds)S. dt

(4)

If we assume the uninfected snails have a constant recruitment rate, A~, then we have from (1), d

-~S = A~ - p~S - phi(P, + P~)S.

(5)

1216

Z, FENC et al.

Combining (2)-(5) and introducing the notation 5i = #h + #p + Oe+ cri,

~ = pbp,

ko# -

k~+l ki

(i = 1, 2),

(6)

we have the six-dimensional system d = Ah,1 -- #hN1 -- C~Pl 4- m2,1N2 - ml,2N1,

7 X* d aix d

ate1

= Ah,2 - #hN2 - aPe + ml,2N1 - ra2,1N2, = /JlrIN1-

61Pl -OlkO l ~ P21~

' \NIJ'

(7)

d

?7P = flarlN2

- 6=P2 - ako,z \ N 2 J

'

d - - S = As - # s S - ~(Pl q- P2)S, dt d Z I = ((P1 + P2)S - (~s + 4 ) I . For mathematical simplicity, we shall make now a simplification t h a t will allow us to carry out an analysis t h a t is not possible in general, namely assume ds = 0. Then, summing up the last two equations in system (7), we have d (S + Z) = As - ~s (S + I). (8) dt We further assume t h a t the total snail population, S + I is at its carrying c a p a c i t y - - s o that S + I is a constant To. Then, we can get rid of equation for S and then have the five-dimensional system d ~-~N1 = Ah,1 -- phN1 - aP1 + rrt2,1N2 - ml,2N1, d ~-~N2 = Ah,2 -- tZhN2 -- aP2 + r~l,2N1 - m2,1N2, d --~jD1 : ~ I T I N 1 - ~ 5 1 P 1 - oJ~0 1 (-jD12~

' \N1]'

-d P2 = fl2rIN2

(9)

a=P=- ako,= \ Np. J '

d --I {(P1 Jr- P2)(To - I) - p s i . dt The existence and the uniqueness of solutions to system (9) can be proved by using standard methods (see, for example, [10]). Moreover, we shall prove t h a t all the variables remain nonnegative and bounded for t > 0 with any nonnegative initial data. Before proving this claim, we need the following lemma. =

LEMMA 1. L e t A, #, a, F, p, and E be nonnegative. I f (xi, x~) is the solution o f the initial value problem d - ~ x l = A - ~ x l - am2,

d

Pxl

pz2

Z 1(0) ~---Xl, 0 > 0,

x2(0) = x2,o > 0, then x l ( t ) > 0 and x2(t) >_ 0 for ali t > O.

(a + e) z l

Schistosomiasis Models

1217

PROOF. First, we note that x~ > 0 when xl > 0 and x2 = 0. This implies xe(t) _> 0 as long as • > 0. By contradiction, we assume that there exists t2 such that xl(te) -- 0, and xl(t) > 0 for

0 < t < te. We now consider two mutually exclusive cases. CASE I. 0 _< me(re) < A / a . From the equation for xl, we know that x~(te) > -#xl(te). By the continuity of me(t), we know there exists to such that x~(t) _> - # x l ( t ) for to _< t _< te. Hence, !

xl _
for to < t < t2. Integrating the above inequality, we have lnxl(t0) - lnxl(t) _< # ( t - - to) <_ #(t2 - to), for to < t < te. As t --~ t~, the left-hand side of the above inequality approaches oo while the right-hand side remains bounded. This means such te does not exist. CASE II.

xe(t2) >_A/a. Consider the set S=_{tc[o,te]lxl(t)=min{X~O),A

pA}} #' Fa

By continuity, we know that S is closed and t2 ~ S. Hence, S has a supremum, to < t2. We have

xl(to) = m i n { X l ? ) h P ~ } P and again, by continuity of xl(t), we have 0 < xl(t) < xl(to) for to < t < te. Note that xz(t) is decreasing for t 6 (t0,t2). If not, from the facts that x~(t) < pA/Fa for to < t < t2, that x[(t) < 0 for t < t2 close enough to t2 and the continuity of the equation for me, we would have the existence of t. such that xe(t,) <_ A / a , but xe(t) > A/a for t. < t < te. However, the fact that x2(t) > A/a for t. < t < t2 and xe(te) > A/a tells us that such t. does not exist, by observing that x~(t) < 0 for t. < t < re. Therefore, it follows that xe(to) > A/a. Next, we shall prove that there exists tl E (to, te) such that xe(tl) = max{xe(to)/2, A/a}. If not, then by continuity, we would have xe(te) > xe(to)/2. We also see from the equation for xe, that

<_ -(a + c) - - , X2

Xl

for all t C (to,te). Integrating from to to te and given that 0 <

x1(t) < xl(to), we have

ln2 _> l n x e ( t 0 ) - lnxe(te) >_ X-z-~ a + e / t l 2 xe(t)dt. Thus, xl(t0)ln2 -> / t l 2 me(t) dt. From the equation for xl and the fact that xl(t) < A/# for t E (t0,t2), we see that x~ _> - a x e . Integrating this relation from to to te, we have Xl(t0) -- Xl(t2) _< a Hence,

:re(t)

( oln )

m1(te)_>xl(to) I

dr. >0.

This contradicts Xl(t2) -- 0. Thus, such h exists. If me(h) > A/a, then we consider tl as a n e w initial point in the above proof. Note that the n e w x2(0) is less than half of the original me(to). Thus, repeating the above procedure a finite number of times, we see that there exists a t* < t2 such that me(t*) -- A / a and xl(t*) > 0. The fact that x2 is decreasing makes me(t2) < A/&. The argument of Case I tells us that is impossible.

Hence, the assumption that such t2 exists is false. The proof is complete. Now, we prove the nonnegativity of solutions of system (9) from nonnegative initial data.

|

1218

Z. FENG et al.

THEOREM 1. If Nl(0) > 0, N2(0) > 0, PI(0) > 0, P2(0) > 0, and I(O) > O, then Nl(t) > O, N2(t) > O, Pl(t) >_ O, P2(t) >_ O, I(t) > 0 for t > O. PROOF. By contradiction, let us assume there exists a ta such that Nl(t3) > 0, N2(t3) > 0, Pl(ta) _> 0, P2(ta) _> 0, I(ta) > 0 but one of these five classes violate the above inequality for t just after t3. We first show that P1 cannot be such class. Indeed, in the equation for P1 in system (9), Pt(ta) = 0, I(t3) _> 0, and Nz(t3) > 0 imply P{(t3) > 0. So, it is easy to see that Pl(t~-) > 0. Thus, we know P1 cannot be one of them. Similar arguments can be applied to P2 and to I. Finally, we claim that N1, ?/2 cannot be such classes either. Without loss of generality, assume that N1 becomes 0 no later than N2. Then, there exists t2 such that N~(t2) = 0 while N~(t) > 0 for t c [0, t2), i = 1,2. We also note that I is bounded above by To, by the nonnegativity of P1 and P2. Now we can generalize Lemma 1 with N1, P1 playing the roles of xl and x2, respectively. Note that replacing A, #, F, p, a + e in Lemma 1 by A h , 1 + rn2,1N2, Ph + ml,2, flirI, 5i, ak0, respectively, in system (9), the only change is that A and F in Lemma 1 become positive functions instead of positive constants. Following the procedure of Lemma 1, we can show that t~he class N2 must violate the positivity inequality earlier than N1, a contradiction. This concludes the proof. | Before proving that the population variables are bounded above, we solve the following system of equations: 0 = Ah,1 - #hN1 + m2,1N2 - ml,2N1, (10) 0 = A~,2 - t~hN2 + rnl,2N1 - m2,1N2. The only solution of (10) is ~N(°)1, N(°)~2j, where

N~o) =

mj,iAh,j + (#h + mj,i)Ah# (#h + mi,j)(#h + mj,i) -- mi,jmj#'

(11)

with i , j = 1,2, and i # j . THEOREM 2. For any nonnegative initial data the solutions of system (9) are bounded for all time. PROOF. In the proof of Theorem 1, we have shown that the class I is bounded above. We will show next that, for i = 1, 2, Ni(t) < N~ °) if Ni(0) < N~ °). If not, then there exists a trajectory that crosses the hyperplane N1 = N} °) or the hyperplane N2 = N (°). Assume, without loss of generality, that the trajectory crosses N1 = N~ °), i.e., there exists a tp such that Nl(tp) = N~ °) and Nl(t +) > N} °) while N2(tp) <_ N (°). Now the equation for N1 says N{(tp) < O, telling us that is impossible. Thus, we see that -/Vi(t) is bounded above by N[ °). Since I and N1 are bounded above, it is easy to see from the equation for P1 in (9) that P1 is bounded above. Similarly, so is P2. | The proofs of existence and uniqueness of solutions, and the well posedness of system (7) are similar to those for system (9). The nonnegativity of S can be proved from Nd S > 0 for S = 0. We also note that

d ( s + I) _<&

+ I).

That gives a upper bound A~/p~ for S + I, and hence, for I and S. 3. A N A L Y S I S

OF

STEADY

STATES

We explore the equilibria of system (9) first. It is easy to see that the parasite-free equilibrium Eo ~ \ rN 1(°) ' N 2(°) , 0, 0, 0) always exists, where N~ °) is given by (11).

Schistosomiasis Models

1219

We introduce now the basic reproductive numbers, 741,742, and 740. Let 74~ = ,/~(~) ~(~) V ' ~MS'~SM'

i = 1, 2,

(12)

where

If we restrict our observation in Village i only, then ~ i plays the same role of the reproductive number 74o in [9] while ~(~) '~MS, ~(i) '~SM correspond to 74MS, 74SM, respectively. In this paper, let ~o = 17421 4- ~22.

(14)

We shall show that 74o plays an important role as threshold for this dynamical system. The stability or instability of the trivial equilibrium determines whether the parasites will be able to establish themselves in an uninfected population. The following result shows that the parasites will go extinct if 74o < 1. PROPOSITION 1. Consider system (9) and assume that P1, P2, and I are twice differentiable with bounded second derivatives. If 74o < 1, then the parasite-free equilibrium Eo is a global attractor, i.e., tlirnoc (N,(t),N2(t),P,(t),P2(t),I(t))

= (N~°),N~°),0, 0,0)

for any positive solution of system (9). In order to prove Proposition 1, we need the /emma below, found in [11]. For a bounded real-valued function f on [0, oo), we define f ~ = liminf f(t),

foo = limsup f(t).

t---+~

t--+ oO

LEMMA 2. Let f : [0, oo) -~ R be bounded and twice differentiable with bounded second derivative. Let t~ --~ oo and f(t,~) converge to f ~ or foo for n ~ oz. Then,

f'(t~) --+ O,

n ~ oo.

PROOF OF PROPOSITION 1. Let 740 < 1. Note that Pl(t), P2(t), I(t) are nonnegative and bounded, and also that 0 < Nl(t) <_ N} °), 0 <_ N2(t) < N (°). Let {sn} be a sequence such that I(sn) I as n - - oo. We can find a subsequenee {sn } of {sn} such that PI,P2 converge for such subsequence. Using Lemma 2 and the equation for I in system (9), we see that

O = ~ ( li+rn Pl(Snk) + liln P2(snk)) ( T o - I ° ° ) - #sI°°, 0 < ( (p~o + p ~ ) To - # , I °°, I °° <

(15)

¢ ( P ~ + P ~ ) To #s

Next, let {tn} be a sequence such that Pl(tn) -+ P ~ as n --+ oo. We can find a subsequence {tn~} of {t,~} such that I, N1 converge for such subsequence. Using again Lemma 2 and the equation for P1 in (9), we see that O = ~ l r (limkk_~I(t~k)

0 <_/31rI°°N~ - 51P~, p~xD < ( / 3 1 T ' ~ iOON~O. -

)

lira Nl(tnk)

-51P~°-akO,llimk__+~Nl(t,~), (16)

Z. FENG et al.

1220

By (15),(16), and the fact that N ~ _< N~ °), we have the following inequality

PF<_

/~lr~ToN} °) ( P ~ + P ~ )

(17)

61#s

Similarly, we have

PF <

~2
-

Summing (17) and (18),

we

(lS)

62#,

have

o _< (PF + P ~ ) ( ~ g - i). Since ~0 < 1 and P ~ , p~o are nommgative, P ~ + P2c~ = 0 is the only possibility to make the above inequality hold. If follows that lim P2(t) = O.

lim Pl(t) = 0, t-'+OO

t--+OO

Hence, limt_+oo I(t) = 0 by (15). Now it is easy to see that lim Nl(t) = N} °),

lim N2(t) = N~ °),

t--+~

t ---+OO

which completes the proof.

|

We have shown that E0 is a global attractor when T~0 < 1. Now, we show that E0 is unstable when T¢0 > 1. PROPOSITION 2. IfT-~ 0 > 1 for (9), then the parasite-free equilibrium Eo is unstable. PaOOF. The characteristic equation of the Jacobian matrix of (9) corresponding to E0 is - [ ( P h + ml,2 + k)(ph + m2,1 + A) -- ml,2m2,1] [a a + a l k 2 + a2A + aa] = 0, where ai = 6i + 62 + #s,

a2 = a , < + (al + a2)~, - eTo~,~N~ °) - ~To92~N~ °), a3 = ai a2>, - {To/31r N~ °) 62 - ~To/32rN~ °) Ji. Since the two roots of the first factor are either negative real or complex conjugates having negative real part, we know by the Routh-Hurwitz criterion that Eo is unstable if a3 < 0, which is equivalent to ~0 > 1. | Proposition 2 suggests that when 7~0 > I the parasites will maintain themselves in the population. In fact, we shall show that result under the assumption that the per capita disease-induced mortality rate of humans, a, is much smaller than the other parameters, which is true for schistosomiasis. PROPOSITION 3. / f R o > 1, then (9) has a unique endemic equilibrium E. = ( N~ , N~ , P~ , P~ , I *) , that is P~, P~, I* > O. Moreover, if a is small, then E . is locally asymptotically stable. PROOF. First, we consider ease a = 0. It is easy to see then that the condition for an endemic equilibrium t3. to exist is ~ o > 1, and we know t77. is unique since

x ; = x~ °),

p; :

~s(~o

- 1)r~

{n~

'

z* -

T°(~°

n2°

- 1)

,

i

=

1, 2.

Schistosomiasis Models To explore the stability of matrix evaluated at E.

E,,

1221

we find the characteristic equation corresponding to the Jacobian

- [ ( . ~ + ~1,2 + a ) ( . ~ + ~2,~ + ~) - ~1,2~:,,1 [ ~ + blX ~ + b~; + b~] = 0,

(19)

where b1=51+52+#~o

2,

ag - nl 2,

b~=,~1,~+<~1~ ~I-~I)+<~.= b3 = 5152ff, (7~ - 1).

Since the two roots of the first factor are either negative real or complex conjugates having negative real part, we know again by the Routh-Hurwitz criterion that E . is locally stable if T~0 > 1 . Let X(t) denote (Nl(t), N2(t), Pl(t), P2(t), I(t)), then we can rewrite system (9) as

x ' = A(~, X).

(20)

We want to show that the nonlinear system A(a, X) = 0 has a unique solution when a is small. We have already pointed out that A(0, X) = 0 has a unique solution E . . We see that the Jacobian matrix of A at (0, E . ) has determinant

--[(]~h -}- 7r~l,2)(]~h + m2,1) -- m1,2m2,1]b3 7d O. By the implicit function theorem in [12], it follows that there exists a unique endemic equilibrium when a is small. Moreover, we also know that such E . is locally asymptotically stable when is small, by the continuity of the eigenvalues of the characteristic equation of E . with respect of a. ! The behavior of system (7) is similar to that of system (9). We know that the parasite-free equilibrium E L = (N}°),N~°),O,O,S(°),O) always exists, where N} 0), N~ 0) are the same as those in system (9), while S (°) = A~/#s. We also define the reproductive numbers, 7~, 7~, and 7Z~)by

~ =,/vU)' V ' " M S ' "~(~)' SM,

(21)

i = 1, 2,

but now

u(2'~ : s (°)

~

,

'~'~ =

, . ~ +d~

"

And

n~ = V / ( ~ ) 2 + ( ~ ) ~ .

(23)

We know that the parasite-free equilibrium E L is a global attractor when T£~ < 1 and it loses its stability when 7¢~ > 1. We also know that the endemic equilibrium E~, is unique when 7Z~ > 1 and a = 0. The implicit function theorem tells us again that E~, is unique when a is small. Though we cannot derive any stability results for El, in this case, the results of ore" simulations suggest E~, is still locally stable when 7Z~)> 1. 4. EXTENSION

TO

n

HUMAN

GROUPS

Since there may be more than two human groups sharing the same water resource, we are going to extend the system from two human groups to n. Using the notation defined above, we

1222

Z. FENG et al.

consider the following systems

j¢i

~p,=

-

j¢i

_

(24)

dt

= ~

P~

S - (.~ + 4)L

i=1

()

and d

j¢i

~dP ~ = 3 i r l N i - 5iPi _

dt

~ £

jT~i

akO'i

\ Ni J

P~ (To - I) - mI,

i=1

where i , j = 1 , . . . , n. The structure of systems (24) and (25) is similar to those of systems (7) and (9), respectively. The existence and uniqueness for these two systems of n human groups are trivial. It is also easy to prove the well posedness of these systems after showing the following 1emma. LEMMA 3. The system

(26)

has a unique positive solution \( N 1(°) ' N,2(°) ' ' ' " , N(°)). PROOF. From elementary linear algebra, it is easy to see that the normal vectors of the above hyperplanes are linearly independent (since all parameters involved are positive, the system is diagonally dominant). Thus, there is a unique solution. If some of those N {°) are nonpositive, say N { ° ) , . . . , N}: ). Then, we see a contradiction by summing the ij th equations for 1 <__j < k, and noting the right-hand side of (26) must be greater than 0. | Let us define the new reproductive numbers of system (25) as

T~,o

I£

7Z2

(27)

i=1

where ~

is given by (12) and (13), while the new reproductive numbers of system (24) are

(2s) where 7~i is given by (21) and (22).

Schistosomiasis Models

1223

Similarly to Proposition 1, it is easy to show that the parasite-free equilibrium is a global attractor if T~,0 < 1 in system (24) or 7~, 0 < 1 in system (25). Since the global stability implies local stability, we have the following corollary by observing the principal n x n submatrix of the (2n + 2) × (2n + 2) Jacobian matrix of system (24) at the boundary equilibrium. COROLLAKY. T h e eigenvalues of the following m a t r i x --~ -- ~ m l , i i~l ml, 2

ml,n

m2,1

...

ran,1

--]_t -- E ?/t2,i i¢2

" ""

?Ttn,2

m2,n

...

(29)

- - # - - ~ mn,i i~n

all have the negative real parts, if # and mi,j are positive real numbers.

We have a result similar to Proposition 2, that is, the parasite-free equilibrium is unstable if 7~n,0 > 1 in system (24) or 7~, o > 1 in system (25). The proof is similar. Summarizing the above results, it is very likely that systems (24) and (25) have structures similar to those of (7) and (9), respectively. Though we cannot prove at present the stability of the endemic equilibria, the simulations in the next section indicate the extension of Proposition 3 to this case should also hold.

5. N U M E R I C A L

SIMULATIONS

In this section, we provide numerical evidence that supports and extends the conclusions of previous sections. In particular, we see the extinction of parasites when the reproductive number is less than unity and we see convergence to an endemic equilibrium when it is greater than unity. First, we employ Newton's method to find the estimations of E . for later comparison with the asymptotic populations reached in the simulations. We find that, to determine E . , we have to solve the following system of equations 0 ---- A3, 0 ( N ~ ) 3 -1- A2,1 (N~') 2 N ~ + . . .

-{- A2, 0 ( N ; ) 2 ÷ . . .

÷ Ao,o,

0 = B3,0 (N~') ~ + B2,1 ( N ; ) 2 N~ + . . . + B2,0 (N~) 2 + " " + B0,0, with explicit--albeit complicated--coefficients A~,js and B~,js. Using Newton's method as described in [13], we approximate N{ and N~ with the following equations 0 = A.,1 -

,

Nf -

+ .

2,1N

-

0 = Ah,2 -- # h N ~ -- o~P~ + m l , 2 N ~ - m2,1N~, 0 = ~ (P{ + P~) (To - I*) - (#s + ds)I*,

while the approximations of P{, P~, and I* are known from Proposition 3. Using a similar method, we obtain an estimation of Et.. Next, we use the backward Euler's method to simulate the dynamics of systems (7) and (9). In order to simulate a fairly realistic situation, we use years as units of time, and we choose the parameters ]_th : 0.014 (human beings average life span = 70 years = 1/#h), AhA ---- 5.6, Ah,2 ----11.2, (without migration, the carrying capacities of Village I and Village 2 are 400 and 800, respectively) #s = 0.2 (snails average life span is 2-8 years), As -- 2000 (without schistosomiasis, the carrying capacity of the snail population is 10,000), k0,1 = 4.31, k0,2 = 4.35, (clumping parameter k is about 0.3, and ko# = (k + 1)/k, as measured from data in Brazil), ds = 0.001, r = 10, ml,2 = 0.0011, m2,t = 0.0009, /~1 = 0.0005, /~2 = 0.001, #p = 0.3 (average life span for a parasite is about three years), a = 0.0000000001, 51,52 = 0.3140000001 or 0.5140000001 (5i = tth ÷ #p + a ÷ (r~, the treatment rate in the first case is 0 while in the second case is 0.2).

1224

Z. FENC et al.

R1<1, R2<1, RO
,

2000 ~

""

0 I , ,

~ ' 50

0

~ ' 100

,-T";'",'i ..................;..............: ' : ; ........................ i ................. ; .......................... ; .................... 150 200 250 300 350 400 450 500 time

R I < I , R2<1, R0
25000

/

jy.~J"

/ 20000

15000

// /

!

10000

',,,

t/`¸¸

5000

0

. . . .

i

. . . .

100

i

. . . .

200

i

. . . .

300

;

. . . .

400

i

. . . .

500

i

. . . .

600

;

. . . .

i

. . . .

i

. . . .

700

800

900

I000

700

800

900

1000

time

R1<1, R2>1, R0>I 80000 70000

// 60000

//

50000 40000 30000 20000 10000

y

0 100

200

300

400

500 lime

Figure 2.

600

Schistosomiasis Models

1225

I~0<1

~*~ II

~p2

i

Ri
1 t~O~

~ooe

0

"100

20~

3;,

,;0

0

~0

~00

tim •

Figure 3.

We use the estimations of E . and E~. to check if (Nl(t), them when TO0 and T¢~, respectively, are greater than 1.

N2(t), Pl(t), P2(t),I(t))

converge to

In Figure 2, we show some results for system (7). We adjust ~ to make the reproductive rate in the range we want. In the first and second graphs, the reproductive numbers 7¢0 are close to 1, so that we can appreciate the importance of this threshold. The first graph in Figure 2 was produced using ~ = 0.0000009, resulting in 7¢1 = 0.55, 7¢2 = 0.82, and T¢0 = 0.99. This, of course, leads to the extinction of all parasites, as clearly indicated by the graph. For the second graph, we used ~ = 0.000000174, giving 7¢1 = 0.49, 7¢~ = 0.94, and 7~0 = 1.06. The corresponding steady-state values for the parasites are P{ = 7,161 and P~ = 26,841, reflected exactly in the simulation. Finally, for the third graph, we used ~ = 0.000001785, that gives 7~1 = 0.77, ~2 = 1.16, and T40 = 1.39. The corresponding steady-state values for the parasites are P{ = 32,093 and P~ = 73,492, again reflected exactly in the simulation. The simulations support the theory about the stability of equilibria for both systems and also suggest that the conjecture about the stability of the endemic equilibrium made for (7) is justified. Moreover, they stress the importance of the transients, since several hundred years are needed in some cases for the system to reach equilibrium and predictions for such length of time are unrealistic. Figure 3 depicts the results of simulations for three human groups. In Figure 3, the parameters used are: Ah,1 = 5.6, Ah,2 = 8.4, Ah, 3 : 11.2, #h ---- 11.2, a = 0.0000000001, rnl,2 = 0.0011, rn2,1 = 0.0009, rnl,a = 0.0012, m3,1 = 0.0008, m2,3 = 0.001, rna,2 = 0.0007, ~1 --'-- 0 . 0 0 0 5 , /~2 = 0 . 0 0 1 , /~3 = 0 . 0 0 1 5 , ~' = 10, 61 = 0.3140000001, 62 = 0.4140000001, 63 = 0.5140000001, k0,1 = 4.31, k0,2 = 4.33, k0,3 = 4.35, #~ = 0.2, As = 2000, and d~ = 0.001. The parameter ~ is 0.00000042 in the upper graph and its 7~0 is 0.9265 while = 0.00000048 in the lower graph and its ~ o is 1.0588. Figure 3 indicates that the discussion in the above section is correct. Though we cannot show (Nl(t), N2(t), Na(t), Pl(t), P2(t), Pa(t), S(t), I(t)) converges to the endemic equilibrium when 7% > 1, choosing the parameter a << 1, we have N[ close to Ni(°) and this is the value approached in the simuiations.

Z. F E N C e t al.

1226

6. C O N T R O L

STRATEGIES

The above results provide a condition, T~0 < 1, for the eradication of the disease. In [8], we see that the treatment rates cq, the snail-to-man transmission parameters ~.i, and the man-to-snail transmission parameter ~ are important in the reduction of ~o- The first, of course, should lead directly to a decrease in the parasite population and, therefore, also in those of miracidia and cercariae. The second one can be decreased by reducing exposure of humans to contaminated waters. The migration parameters rn~j are the new important parameters in this paper. We shall discuss the relationships between these parameters and 7%o. Unless we say otherwise, we adopt the parameters used in the previous numerical simulations.

R0

F i g u r e 4. We choose ~ = 0.0000357, a n d t h e t r e a t m e n t r a t e s ~1 = a2 = 0.

10

6 0

r~l 5

-- - m_21 =0 --

4

m_12=0

\

3 2 1 0

. . . . . '

'

'

t ' '

01

'

I

'

02

'

'

I

03

'

'

'

I

04

'

'

'

I

'

05

'

'

I

'

'

'

06

migrationrate

F i g u r e 5. Z e r o - b o u n d a r y curves.

1'

07

'

'

I

08

'

'

'

I. . . .

0,9

i I

1

Schistosomiasis Models

1227

1.4 1.3 1.2 1.1

R_0 0. 0 (] (

Figure 6.

7.1

0

01

02

0.3

04

05 (72

0.6

Figure 7.

0.7

08

09

1

Z.

1228

/

.......................... __

FENO et al.

.................

~_2

o

o

0

d

F i g u r e 8.

p l

d

d

o

o

c5 d

~

d

d

d

c$ c5 ~

F i g u r e 9.

o

c/

~

d

d

Schistosomiasis Models

1229

Figure 4 shows the relationship between ml,2, rn2,1, and 7~0. Here, we have 7~1 = 3.19 and 7~2 = 6.19. The next figure shows the zero-boundary curves of Figure 4, namely m1,2 = 0 and 7/12,1 ~-~ 0.

We observe in Figures 4 and 5 that T~0 decreases dramatically when ~2,1 is fixed as a small number and ml, 2 increases, and 7~0 increases sharply when mi,2 is fixed as a small number and m2,1 increases. The changes in 7~0 are not so large if the migration rates are large. This shows the migration from the village with the lower 7~i to the one with the higher 7~i helps decrease T~o. On the other hand, the increase of migration from higher 7~i village to lower 7~i village causes an increase of 7~0. However,the effect is less pronounced if the migration rates are large, which in real life they probably are not. In Figures 6 and 7, we show the relationship between or1, a2, and T~0,to understand the relative effect of treatment in each village. We observe the increase of cr2 is much more efficient than that of (71 and the increase of (7] with (72 very small is almost useless. This is just a reflection of the fact that the disease is more prevalent in Village 2. We choose ~ = 0.000001428, so T~0 = 1.39 when the treatment rates (rl, ~r2 are 0. The level curves in Figure 7 show crl as a function of (72 for T~0 equal to 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, and 1.3, respectively. Now, we consider the infection rates/31, and/32 as the control factors. Figure 8 shows 7~0 as a functions of/31 and/32. We choose ~ -- 0.00000357, and the treatment rates (71, ~2 are 0. For each value of 7~0,/31, and/32 have a linear relationship, (23). Figure 9, shows level curves corresponding to T~0 -- 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0, respectively. As we can see in the formula of 7~0, the infection rate of snails ~ -- pbp is proportional to T~0. As we know, there is no practical information on ~. However, without treatment, T~0 ranges from 1 to I00. Substituting the other parameters into the 7~0 formula, we see that ~ ranges from 0.000000357 to 0.0000357. Knowing the range of ~ and the linear relationship with 7~0, we can set up some snail control strategies. The discussion about the extension to n human groups and the resulting reproductive numbers, (27) and (28), tells us that we have to consider every human group surrounding the contaminated water body as a system (instead of a single human group) while investing in schistosomiasis control. Otherwise, the expenses may well be in vain and the effect of treatment rather negligible. On the other hand, if we are in a situation where T~i is extremely small for all i except for T~I, say, then we may reduce 7~0 by controlling T~I only. This control program would save a lot of money.

7. DISCUSSION The mathematical analysis of the two new models in this article derives epidemiological consequences of schistosome distribution among two groups of human hosts, and migrations between human groups. The main contributions of this study are to extend the single human group roodels in [9] to two human habitats, to provide a nearly complete qualitative analysis for the more complex models and to define a threshold that determines the survival of schistosomes. We have shown that, when the reproductive numbers are greater than 1, one can prove the stability of the endemic equilibrimn in system (9) (we obtain it from system (7) by adding an assumption that the snail population is already stabilized). The similarity of these two systems and the results of simulations suggest that the stability of the endemic equilibrium in system (7) is still valid. We also extend the results to n villages, and prove partial results, while the numerical simulations suggest the structures of the models for n human groups are similar to those for the models of two human groups. The nonlinear dependence of T~0 on the treatment rates a~ and on the snail-to-man transmission rates ~ allow us to decide which human population to treat (by increasing its cq) or educate (to reduce its/3i) to design a more efficient treatment program with the available resources.

1230

Z. FENG et al.

We have assumed the release rate of cercariae by one snail is infection-age independent. We make this assumption in order to avoid having a partial differential equation in our system. Control strategies incorporating more empirical data can be discussed as better data becomes available. Future developments of these models can also incorporate seasonal heterogeneity, age-dependent human infection rate, acquired immunity, migration of cerearia driven by flowing waters, and stochastic effects. REFERENCES 1. R.M. Anderson and R.M. May, Regulation and stability of host-parasite population interactions, I, J. Animal Ecology 47, 219-247, (1978), 2. R.M. Anderson and R.M. May, Prevalence of schistosome infections within molluscan populations: Observed patterns and theoretical predictions, Parasitology 79, 63-94, (1979). 3. R.M. Anderson and R.M. May, Helminth infections of humans: Mathematical models, population dynamics, and control, Adv. Parasitol. 24, 1-101, (1985). 4. R.M. Anderson and G.F. Medley, Community control of helminth infection of man and selective chemotherapy, Parasitology 90, 629-660, (1985). 5. P. Hagan and B. Gryseels, Schistosomiasis research and the European community, Tropic. and Geogr. Medic. 46, 259-268, (1994). 6. G. MacDonald, The dynamics of helminth infections with spatial reference to schistosoraes, Trans. Royal Soc. Trop. Med. Hyg. 59, 489-506, (1965). 7. M.E.J. Woolhouse, On the applications of mathematical models of schistosome transmission dynamics h Natural transmission, Acta Tropica 49, 241-270, (i991). 8. M.E.J. Woolhouse, On the applications of mathematical models of schistosome transmission dynamics II: Control, Acta Tropiea 50, 1891-204, (1992). 9. Z. Feng, C. Li and F.A. Milner, Effects of density and age dependence on the transmission dynamics of schistosomes, Math. Biose. 177-178, 271-286, (2002). 10. J.K. Hale, Ordinary Differential Equations, Wiley, (1969). 11. H.R. Thieme, Persistence under relaxed point-dissipativity (with applications to an endemic model), SIAM J. Math. Anal. 24, 407-435, (1993). 12. J.K. Hale and H. Koqak, Dynamics and Bifurcations, Springer-Verlag, (1991). 13. R.L. Burden and J.D. Faires, Numerical Analysis, Fourth Edition~ PWS-KENT Publishing Company, Boston,

MA, (1988).