A Deterministic Model of Schistosomiasis with Spatial Structure

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 for 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 that 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 bilharzia 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.The disease is generally associated with rural poverty.An estimated 170 million people suffered it in sub-Saharan Africa in 2004, and so did a further 30 million in North Africa, Asia, and South America [26].
Schistosomes have to go through an intermediate host (snails in most cases) to complete their life cycle: from eggs, to miracidia, to cercariae, finally to adult flukes.Unlike direct parasites, schistosomes have two stages of reproduction -sexual production in humans and asexual amplification in snails.Mathematical modeling and analysis of schistosomiasis has drawn the attention of many researchers since the first paper by Macdonald in 1965 [15].Thereafter, Anderson and May [2,3], Cohen [6], Nåsell [16], Dobson [7], Adler and Kretzschmar [1], Pugliese [17], and many other researchers built excellent unstructured models and developed a decent understanding of transmission mechanism of schistosomiasis.However, many important biological facts (e.g., age/size-or spatial factors) are overlooked in these models.For example, snails tend to cease production eventually, and they quickly 506 FABIO AUGUSTO MILNER AND RUIJUN ZHAO die after infection [12,20]; children of school age usually exhibit higher prevalence of schistosome infections than other groups [5,24]; there are fewer incidences of infections reported by adventurers who traveled to northern parts of Omo River (faster flow) than by those to southern parts (slower flow) [19].Recently, some age-structured models were proposed [9,24] and optimal control strategies were discussed.However, there is no mathematical model of schistosomiasis with spatial structure.
In this paper, we focus on how the speed of a river affects the transmission dynamics of schistosomiasis by assuming flush-away of only free-swimming miracidia and cercariae, and nothing else.In the literature, lower incidences of schistosomiasis along a river with higher speed occur with fewer snails in the aquatic environment [23].Utzinger et al. believed that Biomphalaria pfeifferi have preference for a certain range of river current velocity, and they claimed a paucity of snails as a reason for low incidences [23].Moreover, spatial microhabitat selection by B. pfeifferi also depends on the water depth [23].On the other hand, when infected snails are prevalent, environmentally triggered downward swimming can quickly bring larvae (planktonic cercariae) to the bed, promoting contact with benthic intermediate hosts [11].Compared with snails, miracidia and cercariae spend most of their life in search for hosts and move along with water flow.Thus, the flush-away of parasite larvae plays an important role in the host-parasite system.Furthermore, we study a control strategy related to treatment and prevention of infected humans and estimate a minimal effort to eradicate the disease.
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 (ODE model, spatially uniform) is studied analytically and numerically; in Section 5, the full model is analyzed and simulations are provided; finally, we present some discussion and conclusions.
2. A mathematical model with spatial structure.Schistosomiasis is found in tropical countries in Africa, the Caribbean, eastern South America, east Asia, and in the Middle East.It is prevalent in villages near rivers or lakes.Since the purpose of our paper is to study the influence of river velocity on the transmission dynamics of schistosomiasis, we assume the aquatic habitat to be a river, which has a geometry of line segment [0, L].Moreover, we assume that the origin of the river (x = 0) is free of parasites, and humans can have an active impact only on a compact region [0, L].
In the following, we make assumptions as in [25].(H0) Human beings and snails grow logistically in the absence of parasites with corresponding carrying capacities.(H1) The rate of parasite-induced human mortality is linearly proportional to the number of parasites a person harbors.(H2) Infected snails do not reproduce and have high potential mortality.(H3) Distribution of parasites among human beings is overdispersed.(H4) Birth rates of hosts are greater than their death rates, b j − µ j0 > 0, j = h, s. (H5) Productivity of cercariae by each infected snail is independent of the dose of miracidia exposed.The biological justification and mathematical methodology can be seen in [2,3,12,14,20,21,22].We further assume that only miracidia and cercariae are spatially dependent.Human beings are active everywhere in the habitat, so that their contribution to the environment can be treated as homogeneous.The same is true for parasites, because they reside inside the human body.We further assume snails are more or less uniformly distributed along the river, because they live in the river bed.
Miracidia and cercariae swim in the river with certain directional preference, such as towards to light [11], and toward intermediate host snails.In this paper, we neglect all these motions and assume that the larvae are flushed away by the river flow.In fact, a cercaria can swim up to 5 meters per day, and a miracidium can swim up to 50 meters per day [13].The speed of larvae swimming is at least a 3order magnitude less than a typical river velocity.It can be ignored.For simplicity, the river velocity is assumed to be a positive constant.
Let H(t), P (t), S(t), I(t) denote the total numbers of human beings, adult parasites, uninfected snails and infected snails at time t, respectively.Let m(x, t), c(x, t) denote the density of free-living miracidia, free-living cercariae at time t and location x, respectively.Let , t)dx denote the total numbers of miracidia and cercariae at time t, respectively.All demographic and epidemiological parameters are listed in Table 1, and most of them are identified in [8,25].
along with homogeneous Dirichlet boundary conditions on the inflow boundary for the two partial differential equations and nonnegative initial distributions Remark 1. System (1) is not self-consistent in that the loss of miracidia and cercariae penetrating into hosts is not reflected in the change of densities of these two larvae.In fact, free-living miracidia and cercariae live at most 48 hours in fresh water, and they lose ability to penetrate into hosts after first several hours when they are active in searching hosts [13].On the other hand, this special treatment allows us to reduce the six-dimensional System (1) to a four-dimensional delay differential equations after explicitly solving the two hyperbolic equations.
Remark 2. Boundary conditions (2) can be modified into some nonnegative functions, in which the model can be applied to study how a contaminated river by schistosomiasis can affect its downstream regions.

3.
Well-posedness.To determine the well-posedness of System (1), we reduce it to an equivalent lower-order dimensional system.Since equations (1.v) and (1.vi) are first-order hyperbolic equations, we can solve them exactly along characteristic lines Since the habitat has finite length L, we can group characteristic lines as above x = vt and below x = vt.The solution of equation ( 4) in the region and the solution of equation ( 4) in the region Thus the solution of equation The total number of miracidia at time t can be obtained by integrating equation ( 7) Similarly, the solution of equation and total number of cercariae at time t is As in [25], we investigate an equivalent system by introducing two new unknowns: mean parasite load per person R = P/H and total number of snail population T = S +I.Substituting equation (8) and equation (10) into System (1) and carrying out some calculation and simplification, we have the following delay differential equations: where M (H(•), R(•)) and C(S(•), T (•)) have the same expressions as in ( 8) and ( 10) except replacing P and I by HR and T − S, respectively.
Lemma 3.1.If H(0), R(0), and S(0) are positive and T (0) > S(0), then solutions of Eqs.(11) It is clear from equation ( 8) that ∀t > 0; if R(s), 0 < s < t is nonnegative, then M (t) is nonnegative.equation (12) implies that if M (t) is nonnegative, then (T − S)(t) is nonnegative.Moreover, from equation (10), if (T − S)(s), 0 < s < t is nonnegative, then C(t) is nonnegative, and equation (11.ii) implies R(t) 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 equation (11.i) gives the positivity of H, equation (11.iii) gives the positivity of S, and then equation (11.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 establish 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 equation (10).Then, R is positive on [0, t] from equation (11.ii) and then equation ( 12) 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 positivity of the solution to System (11).Lemma 3.2.If H(0), R(0), and S(0) are positive and T (0) > S(0), then solutions of equation ( 11) H(t), R(t), S(t) and T (t) are uniformly bounded as long as they exist.

A DETERMINISTIC MODEL OF SCHISTOSOMIASIS WITH SPATIAL STRUCTURE 511
Proof.For convenience, we introduce three new notations: It is evident that equation (11.i) and equation (11.iv) imply that H ≤ L h and T ≤ L s .Moreover, S ≤ L s and I ≤ L s .
For the boundedness of C(t), it suffices to check for t > L/v, From equation (11.ii), we have an estimation The boundedness of R follows immediately Similarly, it suffices to check the boundedness of M (t) for t > L/v, Concerning the existence and uniqueness of solution of the problem, we consult some results by Driver in [4].Consider a system of Volterra functional differential equations Lemma 3.3.(Driver (1962)) Let F(t, x(•) be continuous in t and locally Lipschitz in x and let φ ∈ C([α, t 0 ] → D).Then there exists an h > 0 such that a unique solution x(t) = x(t, t 0 , φ) exists for α ≤ t ≤ t 0 + h.
Rewriting System (11) into a form of equation (17), it is easy to show that F(t, x) is locally Lipschitz in x and is independent of t by observing equation (10), equation (8) and System (11).The local existence and uniqueness of solutions can be obtained by Lemma 3.3.The global existence and uniqueness of solutions is guaranteed by the uniform boundedness.Then, the following theorem has been proved.Theorem 3.4.System (1), or equivalently, System (11) 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.
Remark 3. It is worthy to point out that the domain of φ is {0} for t ≤ L v and [0, L v ] for t > L v when we apply Lemma 3.3 to prove Theorem 3.4.
4. Special case: Model in spatially uniform distribution.At 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 states if nonhomogeneous in-flow boundary conditions are imposed.
Instead of looking at System (1), we study an equivalent system in which the unknowns are H, R, S, T, m and c, It is clear that System (18) admits four parasite-free steady states: where, E DF E represents (H, P, S, I, m, c) = (L h , 0, L s , 0, 0, 0) for System (1).E 0 is the trivial equilibrium, while E 1 and E 2 are semitrivial in that they represent, respectively, the ultimate state of a logistic population of snails and of humans in the absence of parasitism.
Carrying out the similar calculations as in [25], at steady states of System ( 18), H satisfies a cubic equation where Theorem 4.1.Let D, Q and U be defined in equation (20), and consider D and U as functions of treatment rate r.
Remark 4. It can be easily seen from System (18) that if a root H * of equa- ), then System (18) has an interior equilibrium with human population size H * .
As pointed out in [25], it is hopeless to do an analytical stability analysis and so we run some numerical simulations and generate a bifurcation diagram with two parameters L h and r while the rest are fixed.Figure 1  By analyzing the Jacobian matrices of System ( 18) at E 0 , E 1 , E 2 , we can obtain similar results as in [25].
Theorem 4.2.E 0 , E 1 and E 2 are unstable.In region I, there is no interior equilibrium and E DF E is globally asymptotically stable.In region II, there is one interior equilibrium, of which H * corresponds to the smallest root of equation ( 19) in terms of magnitude, and it is globally asymptotically stable.In region III, there is also one globally asymptotically stable interior equilibrium, of which H * corresponds to the largest root.In region IV, there are three interior equilibria, of which the largest and the smallest roots are locally asymptotically stable.In region V, there are also three interior equilibria, of which two largest are unstable so that oscillatory solutions occur.In region VI, there is only one interior equilibrium but unstable, so oscillatory solutions exist.

Define a reproductive number for parasitism
Proof.The Jacobian of System (1) in spatially uniform distribution at the equilibrium It is clear that λ = −(b h − µ h0 ) and λ = −(b s − µ s0 ) are two roots of f (λ) and they are negative from our assumptions.Whether f (λ) has positive roots is determined by the forth order factor  where Following a theorem we owe to Strelitz (1977), the polynomial we have to check the other restrictions to test whether f 1 (λ) is stable or not.Carrying out multiplications and simplifications, it is clear that Moreover notice that Thus, is stable by Strelitz, and so is f (λ).Then, E DF E is locally asymptotical stable.
Remark 5. We conjecture that E DF E is actually a globally asymptotical steady state if R * 0 < 1. Numerical simulations verified that the reproductive number is indeed a threshold for the system.As shown in Figure 4, if R * 0 > 1, the parasitism persists, but if R * 0 < 1, then parasitism disappears.

5.
Model with spatial structure.In Section 3, we established the existence and uniqueness of solution of System (11).In this section we focus on some stability and bifurcation analysis.
From System (1), at steady states, m satisfies  It is not hard to solve m explicitly as Similarly, at steady states, c can be solved as where ξ and η are defined in equation (13).
As shown in previous section, the existence of interior equilibria of System (1) can be established as in Theorem 4.1.Moreover, an analytical stability analysis at these interior equilibria is almost impossible.We turn our attention to the diseasefree equilibria of System (11).Allow us to abuse notations: let E 0 , E 1 , E 1 , E DF E again represent the disease free equilibria of System (11), where where The characteristic polynomial (32) has positive roots at E 0 , E 1 and E 2 under our assumptions; then we prove the conclusion.
Define the reproductive number for System (1) as Proof.It turns out that it is easier to analyze an equivalent delay differential equations of System (11) as follows (34) The linearization of (34) at an equilibrium E DEF gives the following characteristic equation: det The left side of equation ( 36) is a parabola, real part of roots of which are always negative.The right side of equation ( 36) is a decreasing function of λ.It is easy to see that the two functions meet at λ = 0 when R 0 = 1, and meet at a λ > 0 when R 0 > 1. Immediately, this proves that E DF E is unstable when R 0 > 1.
If R 0 < 1, we show that there are no zeros of equation ( 36) with nonnegative real part.Suppose λ * = x 1 + ix 2 to be a complex number with nonnegative real part; i.e., x 1 ≥ 0. It is clear that the modulo of right side of equation (36) evaluated at λ * is greater than or equal to δ(b s + d s ).It is also easy to see that the modulo of the right side of equation (36) evaluated at λ * is less than or equal to βζρb p L h L s ξη.This shows that the modulo of left side of equation ( 36) is always greater than the right side of equation (36) for any λ with λ ≥ 0, when R 0 < 1.It follows immediately that E DF E is locally asymptotically stable when R 0 < 1.
Though we have not proven it yet, E DF E seems to be globally asymptotically stable if R 0 < 1. Remark 6.It is worth noticing that R 0 is a decreasing function with respect to velocity v and treatment rate r.This can be seen by the following: and It is clear that for any v > 0 and L > 0, ξ is a decreasing positive function of v. Noticing that η has same kind of structure as ξ, it is also a decreasing positive function of v.Then, R 0 is a decreasing function of v.It is obvious that it is a decreasing function of r, too.
We run numerical simulations, using the Runge-Kutta Method for O.D.E systems and Crank-Nicolson Characteristic Finite Difference Method for first-order hyperbolic equations.The numerical scheme is of second order.Figure 5 shows profiles of miracidia at nearly steady states for different initial distributions.No matter what initial distribution miracidia start with, the density at steady states is determined by equation (29).Since we assume that human beings, parasites, and snails are spatially independent, the profiles of cercariae at steady states should be similar to those of miracidia.
6. Discussions and conclusions.Endemic schistosomiasis is usually associated with a large human population and relatively slow river flow [19].A large number of human beings provides a big reservoir for destination hosts, and by contaminating the nearby aquatic environment, the residents create a desirable habitat for snails.
We suggest that the carrying capacity for snails should be smaller when the speed of river flow is faster, and logistic growth of hosts in a spatial structured model is necessary.On one hand, fast river flows can flush away snails so that miracidia are inhibited from infecting snails.The asexual reproduction of schistosomiasis could be reduced in this manner.On the other hand, river flows flush away not only snails but also parasite larvae.In fact miracidia and cercariae are very active freeswimming larva they find it easier to follow river flows than do snails.Swimming larvae are drawn to light, so that they strive to stay on the surface of river [11], making it more likely they will be flushed away.The direction of larval swimming is largely vertical [13].This supports our one-dimensional model, in which the spatial motion is due only to the flush-away convection term.
As for control of schistosomiasis, some suggest treating infected humans and preventing further infections in endemic regions; others prefer to curtail snail population either by exterminating the snails or by reducing environmental contamination [18].Both strategies have advantages and disadvantages, and they have worked successfully in controlling or eradicating schistosomiasis in the past [18].In our model, we assume treatment and prevention approaches among humans as the only control method.We proved that the reproductive number decreases as treatment efforts increase, and numerical simulation shows that the mean parasite load per person decreases even with a very small amount of treatment.Based on this, we can estimate a lowest effort to eradicate the disease by making R 0 < 1.Moreover, other control strategies, such as killing snails, can be incorporated into our model with a small modification, and results still hold.Our spatial structured model (1) turns out to be solid.On one hand, the reproductive number for schistosomiasis decreases as the velocity of river flows increases.This can be given as an alternate explanation for why low incidences of schistosomiasis are accompanied by fast river flows.On the other hand, the density of larvae should vary along a river.Figure 5 shows that the density of larvae is low near the origin of a river, but high at the end of the river.
In our model the velocity of river flows is assumed to be a constant.A simple extension of our model is to assume the velocity as a piecewise constant function, which can be caused by different topographies of river beds.To fully understand the real situation, two-dimensional, even three-dimensional model should be proposed, since the distribution of species are different either along a river or across a river.Fully partial differential equations, in which every species are spatially dependent, will be the subject of a further investigation.The main purpose of our model is to provide explanations for biological phenomena, and the model achieved our goal for schistosomiasis.

Figure 1 .
Figure1.Bifurcation Diagram.In region I, there is no interior equilibrium and E DF E is globally asymptotically stable.In region II, there is one interior equilibrium, of which H * corresponds to the smallest root of equation (19) in terms of magnitude, and it is globally asymptotically stable.In region III, there is also one globally asymptotically stable interior equilibrium, of which H * corresponds to the largest root.In region IV, there are three interior equilibria, of which the largest and the smallest roots are locally asymptotically stable.In region V, there are also three interior equilibria, of which two largest are unstable so that oscillatory solutions occur.In region VI, there is only one interior equilibrium but unstable, so oscillatory solutions exist.
together with in-flow boundary condition m(0) = 0.A DETERMINISTIC MODEL OF SCHISTOSOMIASIS WITH SPATIAL STRUCTURE 517

Figure 4 .
Figure 4. Mean parasite load per person for different reproductive number R * 0 .On the left graph, R * 0 = 1.000622; on the right graph, R * 0 = 0.999378, where L h = 800, and r are chosen to match R * 0 .

Figure 5 .
Figure 5. Density plots of miracidia: on the top graph, miracidia start with a sine function; on the bottom graph, miracidia start with a cubic function.Parameters are chosen as L h = 800, v = 1, and the rest are same as in Section 3, except µ m = 1, µ c = 1.