# A Transmission Model of Bilharzia

## A Mathematical Analysis of an Heterogeneous Model

Riveau Gilles<fnref n="1" /> and Sallet Gauthier <fnref n="2" />,<fnref n="3" /> and Tendeng Lena<fnref n="3" />

gilles.riveau@gmail.com ; gauthier.sallet@inria.fr ; lenatendeng@yahoo.fr

**RÉSUMÉ.** On considère un modèle de transmission de la bilharziose prenant en compte les hétérogénéités. Nous calculons le taux de reproduction de base Nous montrons que si $R_0 < 1$, alors l'équilibre sans maladie est globalement asymptotiquement stable. Si $R_0 > 1$, alors il existe un unique équilibre endémique et celui-ci est globalement asymptotiquement stable. Nous considérons ensuite les applications possibles à des données réelles.

**ABSTRACT.** We consider an heterogeneous model of transmission of bilharzia. We compute the basic reproduction ratio $R_0$. We prove that if $R_0 < 1$, then the disease free equilibrium is globally asymptotically stable. If $R_0 > 1$ then there exists an unique endemic equilibrium, which is globally asymptotically stable. We will then consider possible applications to real data

**MOTS-CLÉS** : Bilharziose, modèles de Métapopulations, hétérogénéités spatiales, systèmes forte-ments monotones, stabilité globale et locale.

**KEYWORDS** : Bilharzia, Metapopulation model, Spatial heterogeneity, Strongly monotone systems, Local and global stability.

<footnote n="1">¹ EPLS Saint-Louis</footnote>

<footnote n="2">²INRIA Project MASAIE and IRD UMI UMMISCO and LMAM UMR CNRS 7122</footnote>

<footnote n="3">³INRIA Nancy Grand Est and Gaston Berger University, Saint-Louis Senegal</footnote>
# 1. Introduction

Schistosomiasis, also known as bilharzia after Theodor Bilharz, who first identified the parasite in Egypt in 1851, of all the human parasitic infections, is one of the most widespread. It is second only to malaria in terms of socioeconomic and public health impact in tropical and subtropical areas. An estimated 200 million people in 74 countries are infected with the disease -100 million in Africa alone. School-age children are most likely to become infected with this silent, destructive disease because it is easily contracted while bathing or swimming in water contaminated with the parasite. Children shoulder the majority of schistosomiasis' consequences, especially poor growth and impaired cognitive function. For communities already burdened by poverty and ravaged by scourges such as malaria and HIV/AIDS, schistosomiasis is especially devastating.

Schistosomiasis, or bilharzia, is a parasitic disease caused by trematode flatworms of the genus *Schistosoma*. The cycle of bilharzia is complex.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** *The cycle of Schistosomiasis* (The Carter Centre/Graphic by Al Granberg) [0]</figcaption></figure>

Larval forms of the parasites, *cercariae*, which are released by freshwater snails, penetrate the skin of people in the water. In the body, the larvae develop into adult schistosomes, which live in the blood vessels. The females release eggs, some of which are passed out
of the body in the urine or faeces. Others are trapped in body tissues, causing an immune
reaction.

In urinary schistosomiasis caused by *Schistosoma haematobium*, there is progressive damage to the bladder, ureters and kidneys. In intestinal schistosomiasis, caused by *Schistosoma mansoni*, there is progressive enlargement of the liver and spleen, intestinal damage, and hypertension of the abdominal blood vessels.

Extensive water development has taken place in the north of Senegal over the last decade, resulting in a large increase in the amount of fresh water for irrigation. New freshwater habitats for snail intermediate hosts of schistosomes may be created by the construction of dams and irrigation projects. Greater opportunities for water contact may lead to increases of both urinary and intestinal schistosomiasis, as has been observed in villages along the Senegal river basin. We will work on the data obtained by the research NGO EPLS (Espoir pour la Santé, Saint-Louis, Sénégal, http://www.espoir-sante.org/) in a walo village named Podor, during a campaign in 2004-2005.

Mathematical models are potentially valuable aids to a quantitative understanding of schistosome epidemiology and to the design of control programs. The first model Barbour and others [1, 2] and others have pointed to the potentially strong influence of population and environmental heterogeneities on parasite transmission in endemic communities. This work is a preliminary study of the collaboration with EPLS. A basic theoretical model, inspired from [1], is described that is developed to incorporate the impact heterogeneous transmission rates, and the effects of control measures. We analyze the asymptotic behavior of this model. This work is the first step of the collaboration with EPLS.

## 2. The model

We consider *l* multiple homogeneous populations sites, i.e., patches and *m* water contact sites. We can also differentiate between male and female, simply by considering on the same location two virtual sites : one with males and another with females. This distinction is to account that different age and sex classes used different sites and also that there are different types of activity and the time of day of use. Sex-related patterns vary in relation to behavioural, professional, cultural, and religious factors [10]. Then a human site can be a subpopulation of a precise geographical site.

In the sequel we will use, for these water contact sites, for briefness, the term pond. As said, the human sites can be villages or small clusters of humans. Usually in Senegal, the access to water is not situated in the village. Then a contact water site can either be a location in the bank of the river Senegal where villager take water or be a pond or be a irrigation channel for rice culture ...
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Contact with, and contamination of water. (source PDRSL, ministère de la santé, Sénégal)</figcaption></figure>

We denote, as done in [1], by $w \in \mathbb{R}_+^l$ the vector of mean worm burden of human, each component labelled by its correspondng homogeneous patch and by $y \in \mathbb{R}_+^m$ the vector of the density of infected snails in each water site.

This model is generalization of Macdonald model [3]. MacDonald's model is for one site for human and one location for water. In this Macdonald model the population of humans and snails were assumed to be constant. We assume the analogous hypotheses and denote by $N \in \mathbb{R}_+^m$ the vector of number of snails in each pond.

According to MacDonald [3], the rate of acquisition of egg-laying female schistosomes by the human host population should be proportional to the size of human population, the amount of water contact per person and per day, and the density of infected snails. The human population can usually be supposed to be constant, since the time scale over which it changes is much longer than any of the time scales in the infection process.

Moreover it is assumed, following MacDonald, that the infection of a snail by more than one miracidium does not increase the rate of release of cercariæ. This is supported by experimental studies [2]. This explain why the increase in the worm burden depends only of the density of infected snails and a matrix A which embodies the transmission from snails to man.

If we considerer a human patch i, the variation of the mean number of schistosomes per person is given by

$$ \dot{w}_i = \sum_{j=1}^{m} a_{ij} y_j - \mu_j w_i . $$

Where $y_j$ is the number of patent infected snails in pond j and $a_{ij}$ embodies all the biological process relating to the transmission from snail to human : release of cercariæ, probability that a cercarial encounter leads to an adult parasite, rate of exposure to pond j by a human from site i, the number of snails in pond j. Actually the rate of acquisition should be proportional to the density of infected snails $\frac{y_j}{N_j}$, where $N_j$ is the number of
snails in pond j, we simply incorporate the constant $N_j$ in $a_{ij}$. The parameter $\mu_i$ is the death rate of the female schistosomes in the human host population of site i. This is to take the heterogeneity of each site into account. For example a village, or some subpopulation of a geographical site can have been treated by Praziquantel.

Similarly, if we consider $y_j$ the number of infected snails in pond j, we have

$$ \dot{y}_j = \left(1 - \frac{y_j}{N_j}\right) \sum_{i=1}^{l} b_{ji} w_i - \mu_j y_j $$

Where $b_{ij}$ embodies all the factors relating to transmission from man to snail, from the release of eggs in feces, which is determined by the relation of human of site i to pond j, maturation of sporocysts within the snail. The factor $1 - \frac{y_j}{N_j}$ is simply the proportion of non infected snails. The parameter $\mu_j$ is the death rate of the infected snails in pond j. This depends of the pond, since that either some measures of control can be applied in some ponds (e.g. molluscicide, design of irrigation systems ...), or the ecological situation of the ponds is different.

The model of these $l+m$ equations can be “vectorialized” and written in the following system on $\mathbb{R}_+^{l+m}$

$$ \left\{ \begin{aligned} \dot{w} &= Ay - \operatorname{diag}(\gamma) w \\ \dot{y} &= \operatorname{diag}(\mathbf{1} - \operatorname{diag}(N)^{-1} y) B w - \operatorname{diag}(\mu) y. \end{aligned} \right. \qquad (1) $$

Where $A = (a_{ij})$ is the nonnegative $l \times m$ matrix, which captures the snails to human transmission, $B = (b_{ji})$ is the nonnegative $m \times l$ matrix, which captures the man to snails transmission. The vector $\gamma \in \mathbb{R}_+^m$ is the vector of death rates of female schistosome in each human host population site, and $\mu \in \mathbb{R}_+^m$ is the vector death rates of infected snails in each water site.

We also note by $\mathbf{1}$ the vector $(1, \dots, 1)^T$ of $\mathbb{R}^m$.

The notation $\operatorname{diag}(x)$ represents, if $x \in \mathbb{R}^n$ the $n \times n$ diagonal matrix, whose diagonal terms are given by the components of $x$.

The factor $\operatorname{diag}(\mathbf{1} - \operatorname{diag}(N)^{-1} y) B w$ represents the transmission from man to snails. It is linked to the worm burden $w$ in human sites, and the vector of densities of uninfected snails in the different water sites, which is $\mathbf{1} - \operatorname{diag}(N)^{-1} y$. This link is carried out by the matrix $B$, which embodies the relation of human with water and the biological processes of transmission from man to snails.

To a $n \times n$ matrix $M$ is classically associated a directed digraph $G(M)$ with $n$ vertices [7]: The graph has $n$ vertices. There is an arc from vertex $i$ to vertex $j$ if, and only if the entry of row $i$ and column $j$ of the matrix $M$, $M(i, j) \neq 0$. Conversely to a graph $G$ of order $n$ is associated its adjacency matrix $\operatorname{adj}(G)$.

The symmetry of contact snail-human implies that we have the relation

$$ \operatorname{adj}(G(A)) = \operatorname{adj}(G(B))^T $$
The matrix adj(G(A)) is the adjacency matrix of the digraph associated to A. This is a matrix whose entries are 0 or 1. Then this relation simply says that if there is a connection between a human site and a pond, the humans of the site infect the snails and conversely. The intensity of this connection is given respectively by $A \in \mathbb{R}_+^{l \times n}$ and $B \in \mathbb{R}^{n \times l}$.

The graph of relation between human patches and ponds is then a bipartite graph. Then we can just consider the undirected graph associated to the digraph.

This a consequence of biology, this relation says simply that if a human has contact with water, he can contaminate the pond and also get contaminated.

We will assume, without loss of generality, that the graph of the relations between human and snails, which is given by the following adjacency block matrix

$$ \begin{bmatrix} 0 & \mathrm{adj}(G(A)) \\ \mathrm{adj}(G(B)) & 0 \end{bmatrix} $$

is irreducible :

The nodes of this $l+m$ order graph are the patches and the ponds. There is an arc between a patch and a pond and between a pond and the patch if the humans of the patch visit the pond.

This irreducibility hypothesis is not a loss of generality. Indeed, by symmetry, we can consider the associated undirected graph. For this graph, the connected components coincide with the strong connected components. Then we can just consider the separate connected components of the graph, whose behaviors are independent. Hence our hypothesis is not restrictive.

### 3. Elementary properties of the model

Hereafter we denote $x \le y$ if, for any index $i$, $x_i \le y_i$; $x < y$ if $x \le y$ and $x \ne y$; $x \ll y$ if $x_i < y_i$ for any index $i$;

Since we consider only ponds inhabited with snails (otherwise the pond can be discarded in the model), we have $N \gg 0$.

**Proposition 3.1** *The set*

$$ K = \{(w, y) \in \mathbb{R}_+^{l+m} \mid 0 \le w \le \mathrm{diag}(\gamma)^{-1} N \; ; \; 0 \le y \le N\} $$

is an absorbing compact set on $\mathbb{R}_+^{m+l}$, positively invariant for the system (1).

The proof is straightforward. The model is well posed. The domain $K$ is the biological domain of the model.

**Proposition 3.2**

*System (1) is a strongly monotone system on the open set $\dot{K}$ (open set relatively to $\mathbb{R}_+^{m+l}$)*

$$ \dot{K} = \{(w, y) \mid 0 \le w \ll \mathrm{diag}(\gamma)^{-1} N \; ; \; 0 \le y \ll N\} $$
### Proof

The Jacobian of the system, computed at $(w, y)$, is given by

$$J(w, y) = \begin{bmatrix} -\mathrm{diag}(\gamma) & A \\ \mathrm{diag}(1 - \mathrm{diag}(N)^{-1} y) B & -\mathrm{diag}(\mu) - \mathrm{diag}(B w) \mathrm{diag}(N)^{-1} \end{bmatrix}$$

This matrix is a Metzler matrix (i.e. the off-diagonal terms are nonnegative) and is irreducible, which proves the proposition.

<figure><img src="image_3.png" /></figure>

## 4. The basic reproduction ratio

We use the results and the notations of [4], with the exception that we change the sign $\mathcal{V}$.

The disease free equilibrium (DFE) is $(0, 0)$. The part of the Jacobian, computed at the DFE, coming from the disease transmission is denoted by $F$ and the remaining part is denoted by $V$ (note this is $-V$ in [4]).

$$F = \begin{bmatrix} 0 & A \\ B & 0 \end{bmatrix} \quad V = \begin{bmatrix} -\mathrm{diag}(\gamma) & 0 \\ 0 & -\mathrm{diag}(\mu) \end{bmatrix}$$

The next generation matrix is given by $\mathcal{N} = -FV^{-1}$. The basic reproduction is the spectral radius of the next generation matrix $R_0 = \rho(-FV^{-1}) = \rho(-V^{-1}F)$. Then

$$R_0 = \rho \begin{pmatrix} 0 & \mathrm{diag}(\gamma)^{-1} A \\ \mathrm{diag}(\mu)^{-1} B & 0 \end{pmatrix}$$

### Proposition 4.1

The basic reproduction ratio is the square root of the spectral radius of a $m \times m$ matrix

$$R_0^2 = \rho (\mathrm{diag}(\mu)^{-1} B \mathrm{diag}(\gamma)^{-1} A)$$

This is simply a consequence of the structure of the next generation matrix, and the use of a Cayley reduction.

## 5. Stability Analysis

The following theorem gives the complete stability analysis of system (1)
**Theorem 5.1**

If $\mathcal{R}_0 < 1$ the DFE is globally asymptotically stable. If $\mathcal{R}_0 > 1$ then the DFE is unstable, there exists a unique endemic equilibrium $p \gg 0$ in $\mathbb{K}$, which is globally asymptotically stable.

**Proof**

The proof follows from a theorem of Hirsch [4] improved by Smith [5] combined with the results of Diekmann et al. [6].

From [6] if $\mathcal{R}_0 < 1$ the DFE is locally asymptotically stable and if $\mathcal{R}_0 > 1$ the DFE is unstable. The theorem of Smith says

**Theorem 5.2** Let $F$ be a $C^1$ vector field in $\mathbb{R}^n$, whose flow $\phi$ preserves $\mathbb{R}_+^n$ for $t \ge 0$ and is strongly monotone. Assume that the origin is an equilibrium and that all forward trajectories in $\mathbb{R}_+^n$ are bounded. Suppose that the matrix-valued map $DF: \mathbb{R}^n \to \mathbb{R}^{n \times n}$ satisfies

$$x \gg y \gg 0 \implies DF(x) < DF(y).$$

Then either all trajectories in $\mathbb{R}_+^n$ tend to the origin, or else there is a unique equilibrium $p \gg 0$ which is globally asymptotically stable on $\mathbb{R}_+^n \setminus \{0\}$.

The Theorem is also true for the set $\mathbb{K}$, by adapting easily the Hirsch's proof to $\mathbb{K}$.

The hypotheses of Hirsch-Smith theorem 5.2 are satisfied, for our system, by propositions (3.1) and (3.2). The Jacobian given in the proof of proposition (3.2) satisfies clearly the condition

$$(w_1, y_1) \gg (w_2, y_2) \gg 0 \implies J(w_1, y_1) < J(w_2, y_2).$$

We assume that there exists, in $\mathbb{K}$, an endemic equilibrium $p > 0$. Then $p = (\bar{w}, \bar{y})$ satisfying the equations

$$\left\{ \begin{array}{l} A \bar{y} = \operatorname{diag}(\gamma) \bar{w} \\ \operatorname{diag}(\mathbb{1} - \operatorname{diag}(N)^{-1} \bar{y}) B \bar{w} = \operatorname{diag}(\mu) \bar{y} \end{array} \right. \qquad (2)$$

Then

$$\operatorname{diag}(\mathbb{1} - \operatorname{diag}(N)^{-1} \bar{y}) B \operatorname{diag}(\gamma)^{-1} A \bar{y} = \operatorname{diag}(\mu) \bar{y}$$

Equivalently, using the matrix $R = \operatorname{diag}(\mu)^{-1} B \operatorname{diag}(\gamma)^{-1} A$, we have, using the commutativity of diagonal matrices.

$$\bar{y} = \operatorname{diag}(\mathbb{1} - \operatorname{diag}(N)^{-1} \bar{y}) R \bar{y} < Ry$$

The strict inequality comes from $\operatorname{diag}(N)^{-1} \bar{y} \ll \mathbb{1}$.

It is a simple consequence of irreducibility, that $R$ is a nonnegative irreducible matrix. Then we have $y < Ry$, from [7] Theorem 3.31, we deduce $\mathcal{R}_0^2 = \rho(R) > 1$. The
DFE is unstable. By a theorem of Varga, this implies that the Jacobian, computed at the DFE is unstable, with the stability modulus positive. The Jacobian at the DFE is given by $J = F + V$. If we denote the stability radius of $J$ by $s(J)$, by Perron-Frobenius theorem we deduce that there exists $v \gg 0$ such that

$$Jv = s(J)v$$

Then when $R_0 \le 1$ there is no endemic equilibrium. Hence by Smith result the DFE is globally asymptotically stable.

When $R_0 > 1$ then the DFE is unstable, and there exists an unique endemic equilibrium $p = (\bar{w}, \bar{y}) \gg 0$ in $\bar{K}$, which attracts all the trajectories.

Moreover, since $K$ is absorbing, this proves the global stability of $p$.

<figure><img src="image_4.png" /></figure>

## 6. Stability of the DFE when $R_0 = 1$

During the CARI conference the question of “what happens when $R_0 = 1$” has been posed. We address in this section this issue

### Theorem 6.1

If the network is strongly connected, the DFE is globally asymptotically stable on the domain $\mathbb{R}^{mp+n} \times [0, 1]^n$, when $R_0 = 1$.

### Proof

We use the equivalence between $R_0 = 1$ and $s(J) = 0$, where $s(J)$ denotes the spectral bound of the matrix $J$, i.e., the largest real part of the eigenvalues of $J$ (see for example [8, 4]. But $J$ is an irreducible Metzler matrix.

Hence there exists a positive vector $v = (\alpha, \beta) \gg 0$ such that $J^T v = s(J)v = 0$ [9].

We have then the following relations for $\alpha \gg 0$ and $\beta \gg 0$

$$ \left\{ \begin{array}{l} A^T \alpha - \operatorname{diag}(\mu) \beta = 0 \\ -\operatorname{diag}(\gamma) \alpha + B^T \beta = 0. \end{array} \right. \qquad (3) $$

We now consider, on the attracting compact $K$ defined in proposition 3.1, the following radially unbounded Lyapunov function

$$V(w, y) = \langle w | \alpha \rangle + \langle y | \beta \rangle = \left\langle \begin{bmatrix} w \\ y \end{bmatrix} \middle| \begin{bmatrix} \alpha \\ \beta \end{bmatrix} \right\rangle$$

We write system (1) “vectorially”

$$\begin{bmatrix} \dot{w} \\ \dot{y} \end{bmatrix} = M(w, y) \begin{bmatrix} w \\ y \end{bmatrix}.$$

Where the matrix $M(w, y)$ is defined by
$$
M(w, y) = \begin{bmatrix} -\mathrm{diag}(\gamma) & A \\ \mathrm{diag}(1 - \mathrm{diag}(N)^{-1} y) B & -\mathrm{diag}(\mu) \end{bmatrix}.
$$

We observe that this matrix is related to the Jacobian $M(w, y) = J(0, y)$.

The derivative $\dot{V}$ of $V$ along the trajectories of system (1) is given by

$$
\dot{V} = \left\langle M(w, y) \begin{bmatrix} w \\ y \end{bmatrix} \middle| \begin{bmatrix} \alpha \\ \beta \end{bmatrix} \right\rangle = \left\langle \begin{bmatrix} w \\ y \end{bmatrix} \middle| M(w, y)^T \begin{bmatrix} \alpha \\ \beta \end{bmatrix} \right\rangle
$$

Then

$$
\dot{V} = \langle -\operatorname{diag}(\gamma) \alpha + B^T \beta - B^T \operatorname{diag}(N)^{-1} \operatorname{diag}(y) \beta | w \rangle + \langle A^T \alpha - \operatorname{diag}(\mu) \beta | y \rangle
$$

Using the relations (3)

$$
\dot{V} = -\langle B^T \operatorname{diag}(N)^{-1} \operatorname{diag}(y) \beta | w \rangle = -\langle \operatorname{diag}(N)^{-1} \operatorname{diag}(y) \beta | Bw \rangle \le 0
$$

This proves the stability of the DFE. To prove the global stability we will use Lasalle's
invariance principle on $K$.

We consider the largest invariant set $\mathcal{L}$ contained in the set $E$

$$
E = \{ (w, y) \in K \mid \dot{V} = \langle \operatorname{diag}(N)^{-1} \operatorname{diag}(y) \beta | Bw \rangle = 0 \}
$$

We assume first that $y \gg 0$, with $(w, y) \in \mathcal{L}$. Since $\beta \gg 0$, then we deduce immedia-
tely that $\operatorname{diag}(N)^{-1} \operatorname{diag}(y) \beta \gg 0$.

Then to be in $\mathcal{L}$, with $y \gg 0$, implies $Bw = 0$. The system for $y$ reduces to

$$
\dot{y} = -\operatorname{diag}(\mu)y,
$$

which is an asymptotically stable linear system. Any trajectory followed backward
from an initial point (w, y) with y ≫ 0 will leaves any compact K. Then no point (w, y),
with y ≫ 0, can be in L.

Then we assume $y \gg 0$. Let $j$ an index such that $y_j = 0$. Since $y \in \mathcal{L}$ and $y_j = 0$ we
have, by invariance,

$$
\dot{y}_j = \sum_{i=1}^{l} b_{j,i} w_i = 0
$$

This implies that $w_i = 0$ for any $b_{j,i} \neq 0$. In other words this means that the indi-
viduals from the patch of index $i$ (ranging from 1 to $l$), which are excreting in patch $j$
(because $b_{j,i} \neq 0$), necessarily are not infected (i.e. $w_i = 0$) This is intuitively evident.

Now for a $w_i = 0$, in $\mathcal{L}$, we have

$$
\dot{w}_i = \sum_{k=1}^{m} a_{j,k} y_k
$$
We have, again by invariance, that $y_k = 0$ for any $a_{j,k} \neq 0$. In other words, all the ponds, which are visited by the individuals of patch whose index is $j$, must be without infected snails. This is again intuitively evident.

By abuse of language, we identify the node of the graph and the value $y_j$ or $w_i$. Now by irreducibility, starting from a $y_j = 0$ we can access, following a directed path on the graph, any $y_k$. Then we have proved that for all index $k$, we have $y_k = 0$, i.e., $y = 0$. Again, by irreducibility any $w_i$ is accessible from a $y_j$. From what is preceding, this proves $w = 0$.

Then $\mathcal{L} = \{0\}$. By Lasalle's invariance principle the DFE is globally asymptotically stable on $K$. Since $K$ is an attracting set this proves the global stability on the biological domain.

■

## 7. Numerical simulation

We conduct some numerical simulations, with 8 patches and 3 ponds, using SCILAB. The connectivity matrix are

A=

<table><tr><td>0.0001389</td><td>0.0002778</td><td>0.</td></tr><tr><td>0.0016667</td><td>0.0020833</td><td>0.</td></tr><tr><td>0.0018056</td><td>0.0022222</td><td>0.</td></tr><tr><td>0.0013889</td><td>0.0013889</td><td>0.</td></tr><tr><td>0.</td><td>0.0002778</td><td>0.0000694</td></tr><tr><td>0.</td><td>0.0020833</td><td>0.0008333</td></tr><tr><td>0.</td><td>0.0022222</td><td>0.0013889</td></tr><tr><td>0.</td><td>0.0013889</td><td>0.0011111</td></tr></table>

and

B =

<table><tr><td>0.0076</td><td>0.0912</td><td>0.0988</td><td>0.076</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td></tr><tr><td>0.0152</td><td>0.114</td><td>0.1216</td><td>0.076</td><td>0.0152</td><td>0.114</td><td>0.1216</td><td>0.076</td></tr><tr><td>0.</td><td>0.</td><td>0.</td><td>0.</td><td>0.0038</td><td>0.0456</td><td>0.076</td><td>0.0608</td></tr></table>

For exposition we assume that the mortality for worms and snails is always the same, equal to 0.0092593 for the worms and for the snails to 0.038. The ponds 2 and 3 have been cleaned from their infected snails.

We have two groups of 4 villages. The first group uses ponds 1 and 2 and the second group uses ponds 2 and 3. The basic reproduction ratio is obtained readily by the SCILAB command

R02=max(real(spec inv(diagMU)*B*inv(diagGAMMA)*A)))
ans =
4.4419854
We obtain the following curves

<figure><img src="image_5.png" /><figcaption>Figure 3. Simulation of 8 sites and 3 ponds</figcaption></figure>

Now we treat with PZQ the 3 villages during one month, increasing by 10 the death rates of worms. The basic reproduction ratio is then divided by 10 and $R_0^2 = 0.44$. We obtain the following curves on 150 days, starting from the same initial point as the example of the preceding curves. We represent only patches 1-4, since the curves are analogous for the patches 5-8.

<figure><img src="image_6.png" /><figcaption>Figure 4. Simulation of a one month treatment</figcaption></figure>

## 8. Discussion

The chemotherapy of population acts on diag(γ), a biological control will acts on diag(μ). With the value of $R_0$

$$R_0^2 = \rho (\text{diag}(\mu)^{-1} B \text{diag}(\gamma)^{-1} A)$$
and the properties of nonnegative matrices, augmenting the vector $\gamma$ will reduce $R_0$. But the relation is far from linear, if we differentiate the intervention on the patches. We have now to explore numerically the most effective way to reduce $R_0$.

This suppose that the intervention, in this case chemotherapy of some population, is continuous. In reality the intervention are discontinuous. The model must be extended to "pulse" chemotherapy. This will be addressed elsewhere.

We thank the Carter Center and The ministère de la santé du Sénégal for granting per-mission to include the graphics in figure 1 and 2.

## 9. Bibliographie

[0] http://www.cartercenter.org/health/schistosomiasis/index.html

[1] BARBOUR, A. D., « Macdonald's model and the transmission of bilharzia. », *Trans R Soc Trop Med Hyg*, vol. 78, n° 1, pp 6–15, 1978.

[2] WOOLHOUSE, M E, « Mathematical models of transmission dynamics and control of schistosomiasis. », *Am J Trop Med Hyg*, vol. 55, n° 5 Suppl, pp 144–148, 1996.

[3] MACDONALD, G, « The dynamics of helminth infections, with special reference to schistosomes. », *Trans R Soc Trop Med Hyg*, vol. 59, n° 5 Suppl, pp 489–506, 1965.

[4] VAN DEN DRIESSCHE, P. AND WATMOUGH, J., « reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission », *Math. Biosci.*, vol. 180, pp 29–48, 2002.

[4] HIRSCH, M.W., « The dynamical system approach to differential equations », *Bull. AMS*, vol. 11, n° 1, pp 1–64, 1984.

[5] SMITH, H. L., « Cooperative systems of differential equations with concave nonlinearities », *Nonlinear Anal.*, vol. 10, n° 10, pp 1037–1052, 1986.

[6] DIEKMANN, O. AND HEESTERBEEK, J. A. P. AND METZ, J. A. J., « On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations », , vol. 28, n° 4, pp 365–382, 1990.

[7] BERMAN, A. AND PLEMMONS, R. J., « Nonnegative matrices in the mathematical sciences », , SIAM, Classics in Applied Mathematics 1994.

[8] KAMGANG, J.C. AND SALLET, G., « Computation of threshold conditions for epidemiological models and global stability of the disease free equilibrium. », *Math. Biosci.*, vol. 213, n° 1, pp 1–12, 2008.

[9] SMITH, H.L., « Monotone dynamical systems. An introduction to the theory of competitive and cooperative systems. », *AMS*, 1995.

[10] GRYSEELS AND POLMAN, K. AND CLERINX, J. AND KESTENS, L., « Human schistosomiasis », *Lancet*, vol. 368, n° 9541, pp 1106–1118, 2006.