The role of disease cycles in the endemicity of infectious diseases

TL;DR

This paper presents a simple epidemiological model demonstrating how reservoir cycles influence disease endemicity, highlighting that focusing solely on vectors or reservoirs may be insufficient for disease control.

Contribution

The study introduces a network-based model with a bifurcation parameter to explain disease endemicity driven by reservoir cycles, even when individual reservoir cycles are weak.

Findings

Reservoir cycles with small reproductive numbers can sustain endemicity.

Endemic states can occur without strong reservoir cycles when considering the full reservoir network.

Control strategies should target both vectors and reservoirs, not just one.

Abstract

Vector-borne diseases with reservoir cycles are complex to understand because new infections come from contacts of the vector with humans and different reservoirs. In this scenario, the basic reproductive number $\mathcal{R}^h_0$ of the system where the reservoirs are not included could turn out to be less than one, yet, an endemic equilibrium be observed. Indeed, when the reservoirs are taken back into account, the basic reproductive number $\mathcal{R}_0^r$, of only vectors and reservoirs, explains the endemic state. Furthermore, reservoirs cycles with a small basic reproductive number could contribute to reach an endemic state in the human cycle. Therefore, when controlling for the spread of a disease, it could not be enough to focus on specific reservoir cycles or only on the vector. In this work, we created a simple epidemiological model with a network of reservoirs where $\mathcal{R}_0^r$ is a bifurcation parameter of the system, explaining disease endemicity in the absence of a strong reservoir cycle. This simple model may help to explain transmission dynamics of diseases such as Chagas, Leishmaniasis and Dengue.

Tables

Table 1. Table 1: Parameters of the model described by equations ( 2 ).

Parameter	meaning	Units
${\beta }_{vh}$	Number of human infections caused by	$\left[H\right]/\left(\left[time\right]\ast \left[V\right]\right)$
	one infectious vector per unir of time
${\beta }_{hv}$	Number of vector infections caused by	$\left[V\right]/\left(\left[time\right]\ast \left[H\right]\right)$
	one infectious human per unir of time
${\beta }_{v{r}_i}$	Number of infections of reservoir $R_i$ caused by	$\left[R_i\right]/\left(\left[time\right]\ast \left[V\right]\right)$
	one infectious vector per unir of time
${\beta }_{r_{i}v}$	Number of vector infections caused by	$\left[V\right]/\left(\left[time\right]\ast \left[R_i\right]\right)$
	one infectious reservoir $R_i$ per unir of time
${\gamma }_{ij}$	Number of infections of reservoir $R_j$ caused by	$\left[R_j\right]/\left(\left[R_i\right]\ast \left[time\right]\right)$
	one infectious reservoir $R_i$ per unir of time
${\delta }_h$	Mortality rate of humans	$1/\left[time\right]$
${\delta }_v$	Mortality rate of vectors	$1/\left[time\right]$
${\delta }_{r_i}$	Mortality rate of reservoirs $R_i$	$1/\left[time\right]$

Table 2. Table 2: Umbral values defined from the systems formed by the state variables of the second column and the interpretations of the next generation matrix given by the third column.

Value	System form by	New infections
${ℛ}_0$	$H,V,R$	$H$
${ℛ}_0^r$	$V,R$	$V$
${ℛ}_0^h$	$H,V$	$H$

Table 3. Table 3: Umbral values defined from the systems formed by the state variables of the second column and the interpretations of the next generation matrix given by the third column

Value	System form by	New infections
${ℛ}_0$	$H,V,R_i$, $i=1,\dots ,k$	$H$
${ℛ}_0^r$	$V,R$, $i=1,\dots ,k$	$V$
${ℛ}_0^h$	$H,V$	$H$
${ℛ}_0^{r_i}$	$V,R_i$	$V$

Table 4. Table 4: Parameters of the Chagas model in [ 2 ] .

Parameter	Units	Estimate
${\beta }_{hv}^{\prime }$	Fraction of vectors infected by	$1/\left(year\ast \left[H\right]\right)$	0-1
	one infectious human per year
${\beta }_{vh}^{\prime }$	Fraction of humans infected by	$1/\left(year\ast \left[V\right]\right)$	${\beta }_{hv}^{\prime }$/100
	one infectious vector per year
${\beta }_{r_{1}v}^{\prime }$	Fraction of vectors infected by	$1/\left(year\ast \left[R_1\right]\right)$	2 ${\beta }_{hv}^{\prime }$
	one infectious $R_1$ host per year
${\beta }_{v{r}_1}^{\prime }$	Fraction of $R_1$ hosts infected by	$1/\left(year\ast \left[V\right]\right)$	${\beta }_{hv}^{\prime }$/10
	one infectious vector per year
${\beta }_{r_{2}v}^{\prime }$	Fraction of vectors infected by	$1/\left(year\ast \left[R_2\right]\right)$	1 ${\beta }_{hv}^{\prime }$
	one infectious $R_2$ host per year
${\beta }_{v{r}_2}^{\prime }$	Fraction of $R_2$ hosts infected by	$1/\left(year\ast \left[V\right]\right)$	${\beta }_{hv}^{\prime }$/5
	one infectious vector per year
$N_{r_1}$	Number of $R_1$ individuals	$\left[R_1\right]$	0.0005 $N_v$
$N_{r_2}$	Number of $R_2$ individuals	$\left[R_2\right]$	0.001 $N_v$
$N_h$	Number of humans	$\left[H\right]$	0.001 $N_v$
${\delta }_v$	Mortality rate of vectors	1/year	1
${\delta }_{r_1}$	Mortality rate of $R_1$ hosts	1/year	0.5
${\delta }_{r_2}$	Mortality rate of $R_2$ hosts	1/year	0.3
${\delta }_h$	Mortality rate of humans	1/year	0.015
${\beta }_{vh}$	Number of human infections caused by	$\left[H\right]/\left(years\ast \left[V\right]\right)$	$N_h{\beta }_{vh}^{\prime }$
	one infectious vector per unir of time
${\beta }_{hv}$	Number of vector infections caused by	$\left[V\right]/\left(years\ast \left[H\right]\right)$	$N_v{\beta }_{hv}^{\prime }$
	one infectious human per unir of time
${\beta }_{v{r}_i}$	Number of infections of reservoir $R_i$ caused by	$\left[R_i\right]/\left(years\ast \left[V\right]\right)$	$N_{r_i}{\beta }_{v{r}_i}^{\prime }$
	one infectious vector per unir of time
${\beta }_{r_{i}v}$	Number of vector infections caused by	$\left[V\right]/\left(years\ast \left[R_i\right]\right)$	$N_v{\beta }_{r_{i}v}^{\prime }$
	one infectious reservoir $R_i$ per unir of time

Table 5. Table 5: Basic reproductive numbers of the system of equations ( 1 ). We define μ h := ν h δ h assign subscript 𝜇 ℎ subscript 𝜈 ℎ subscript 𝛿 ℎ \mu_{h}:=\frac{\nu_{h}}{\delta_{h}} and μ v ​ h := β h ​ v δ h ​ β v ​ h δ v assign subscript 𝜇 𝑣 ℎ subscript 𝛽 ℎ 𝑣 subscript 𝛿 ℎ subscript 𝛽 𝑣 ℎ subscript 𝛿 𝑣 \mu_{vh}:=\sqrt{\frac{\beta_{hv}}{\delta_{h}}\frac{\beta_{vh}}{\delta_{v}}} .

Interpretation	$F$	$V$	$G=F{V}^{-1}$	$\rho \left(F{V}^{-1}\right)$
(New infections)
$H$				${ℛ}_0^2={\mu }_h+{\mu }_{vh}^2$
$V$				${ℛ}_0^3={\mu }_{vh}^2/\left(1-{\mu }_h\right)$
$H,V$				${ℛ}_0^1={\lambda }^{\ast }$, where
				${\left({\lambda }^{\ast }\right)}^2-\left({\mu }_h\right){\lambda }^{\ast }-{\mu }_{vh}^2=0$.

Equations

$$\{\begin{array}{l}\dot{I_h}={\beta }_{vh}{I}_v\frac{S_h}{N_h}+{\nu }_h{I}_h\frac{S_h}{N_h}-{\delta }_h{I}_h\\ \dot{I_v}={\beta }_{hv}{I}_h\frac{S_v}{N_v}-{\delta }_v{I}_v\end{array}$$

$$\{\begin{array}{l}\dot{I_h}={\beta }_{vh}{I}_v\frac{S_h}{N_h}-{\delta }_h{I}_h\\ \dot{I_{r_i}}={\beta }_{v{r}_i}{I}_v\frac{S_{r_i}}{N_{r_i}}+{\nu }_i{I}_{r_i}\frac{S_{r_i}}{N_{r_i}}+{\sum }_{j\ne i}{\gamma }_{ji}{I}_{r_j}\frac{S_{r_i}}{N_{r_i}}-{\delta }_{r_i}{I}_{r_i},i=1,\dots ,k\\ \dot{I_v}={\beta }_{hv}{I}_h\frac{S_v}{N_v}+{\sum }_{i=1}^k{\beta }_{r_{i}v}{I}_{r_i}\frac{S_v}{N_v}-{\delta }_v{I}_v\end{array}$$

$$\Omega :=\{\sigma :\{i_1\dots i_k\}\to \{i_1,\dots ,i_k\}:\sigma \text{ is a cycle},\{i_1,\dots ,i_k\}\subseteq \{1,\dots ,n\}\}.$$

$$li{m}_{r\to \infty }\mu (A):=\rho (A^{[r]}{)}^{1/r}=ma{x}_{\sigma \in {\Omega }_n}{\mu }_{\sigma }(A)$$

$$\mu (A)\le \rho (A)\le \rho (S(A))\mu (A).$$

$${\mu }_{\ast }(G)\le {\mathcal{R}}_0\le \rho (S(G)){\mu }_{\ast }(G)$$

$$F=\left(\begin{array}{cccccc}0& {\beta }_{vh}& 0& 0& \dots & 0\\ {\beta }_{hv}& {\nu }_v& {\beta }_{r_{1}v}& {\beta }_{r_{2}v}& \dots & {\beta }_{r_{k}v}\\ 0& {\beta }_{v{r}_1}& {\nu }_{r_1}& {\gamma }_{21}& \dots & {\gamma }_{k1}\\ 0& {\beta }_{v{r}_2}& {\gamma }_{12}& {\nu }_{r_2}& \dots & {\gamma }_{k2}\\ ⋮& \ddots & \ddots & \ddots & \dots & {\gamma }_{kk-1}\\ 0& {\beta }_{v{r}_k}& {\gamma }_{1k}& {\gamma }_{2k}& \dots & {\nu }_{r_k}\end{array}\right),V=\left(\begin{array}{cccccc}{\delta }_h& 0& 0& 0& \dots & 0\\ 0& {\delta }_v& 0& 0& 0& 0\\ 0& 0& {\delta }_{r_1}& 0& 0& 0\\ 0& 0& 0& {\delta }_{r_2}& \dots & 0\\ ⋮& \ddots & \ddots & \ddots & \dots & 0\\ 0& 0& 0& 0& \dots & {\delta }_{r_k}\end{array}\right)$$

$$G=F{V}^{-1}=\left(\begin{array}{cccccc}0& \frac{{\beta }_{vh}}{{\delta }_v}& 0& 0& \dots & 0\\ \frac{{\beta }_{hv}}{{\delta }_h}& \frac{{\nu }_v}{{\delta }_v}& \frac{{\beta }_{r_{1}v}}{{\delta }_{r_1}}& \frac{{\beta }_{r_{2}v}}{{\delta }_{r_2}}& \dots & \frac{{\beta }_{r_{k}v}}{{\delta }_{r_k}}\\ 0& \frac{{\beta }_{v{r}_1}}{{\delta }_v}& \frac{{\nu }_{r_1}}{{\delta }_{r_1}}& \frac{{\gamma }_{21}}{{\delta }_{r_2}}& \dots & \frac{{\gamma }_{k1}}{{\delta }_{r_k}}\\ 0& \frac{{\beta }_{v{r}_2}}{{\delta }_v}& \frac{{\gamma }_{12}}{{\delta }_{r_1}}& \frac{{\nu }_{r_2}}{{\delta }_{r_2}}& \dots & \frac{{\gamma }_{k2}}{{\delta }_{r_k}}\\ ⋮& \ddots & \ddots & \ddots & \dots & \frac{{\gamma }_{kk-1}}{{\delta }_{r_k}}\\ 0& \frac{{\beta }_{v{r}_k}}{{\delta }_v}& \frac{{\gamma }_{1k}}{{\delta }_{r_1}}& \frac{{\gamma }_{2k}}{{\delta }_{r_2}}& \dots & \frac{{\nu }_{r_k}}{{\delta }_{r_k}}\end{array}\right)$$

$$\lambda (({\nu }_v/{\delta }_v-\lambda )-\sum _{i=1}^n\frac{{\mu }_{v{r}_i}^2}{{\nu }_{r_i}/{\delta }_{r_i}-\lambda })=-{\mu }_{vh}^2$$

$${\mathcal{R}}_0=\rho (F{V}^{-1})=\sqrt{\frac{{\beta }_{vh}}{{\delta }_v}\frac{{\beta }_{hv}}{{\delta }_h}+\sum _{i=1}^k\frac{{\beta }_{v{r}_i}}{{\delta }_v}\frac{{\beta }_{r_{i}v}}{{\delta }_{r_i}}}=\sqrt{{\mu }_{vh}^2+{\mu }_{v{r}_1}^2+\dots +{\mu }_{v{r}_k}^2}$$

$${\mu }_{\ast }^2\le {\mu }_{vh}^2+{\mu }_{v{r}_1}^2+\dots +{\mu }_{v{r}_k}^2\le (k+1){\mu }_{\ast }^2$$

$${\mu }_{\ast }=max\{{\mu }_{vh},{\mu }_{v{r}_1},\dots ,{\mu }_{v{r}_k},{\nu }_v/{\delta }_v,{\nu }_{r_1}/{\delta }_{r_1},\dots ,{\nu }_{r_k}/{\delta }_{r_k}\}.$$

$${\mu }_{\ast }\le {\mathcal{R}}_0\le \rho (S){\mu }_{\ast }=(\sqrt{k+1}+1){\mu }_{\ast }$$

$${\mathcal{R}}_0=\frac{{\mathcal{R}}_0^h}{1-{\mathcal{R}}_0^r}$$

$${\mathcal{R}}_0^r=\sum _{j=1}^k{\mathcal{R}}_0^{r_j}$$

$$N_{r_2}>-\frac{{\mu }_{v{r}_1}^{\prime }}{{\mu }_{v{r}_2}^{\prime }}{N}_{r_1}+(\frac{1}{N_v{\mu }_{v{r}_2}^{\prime }}-\frac{N_h{\mu }_{vh}^{\prime }}{{\mu }_{v{r}_2}^{\prime }})=-\frac{3}{5}{N}_{r_1}+(\frac{3}{2N_v{\beta }_{hv}^2}-N_h)$$

$$\frac{d{x}_i}{dt}=f_i(\bar{x},\bar{y})-v_i(\bar{x},\bar{y}),i=1,\dots ,n.$$

$${\mathcal{R}}_0^h=\frac{{\beta }_{hv}{\beta }_{vh}}{{\delta }_h{\delta }_v}$$

$${\mathcal{R}}_0^{r_i}=\frac{{\beta }_{r_{i}v}{\beta }_{v{r}_i}}{({\delta }_{r_i}-{\nu }_i{N}_{r_i}){\delta }_v}$$

$$F=\left(\begin{array}{cccccc}0& {\beta }_{vh}& 0& 0& \dots & 0\\ 0& 0& 0& 0& \dots & 0\\ ⋮& & & \ddots & \ddots & ⋮\\ 0& 0& 0& 0& \dots & 0\\ ⋮& & & \ddots & \ddots & ⋮\\ 0& 0& 0& 0& \dots & 0\end{array}\right),V=\left(\begin{array}{cccccc}{\delta }_h& 0& 0& 0& \dots & 0\\ -{\beta }_{hv}& {\delta }_v& -{\beta }_{r_{1}v}& -{\beta }_{r_{2}v}& \dots & -{\beta }_{r_{k}v}\\ 0& -{\beta }_{v{r}_1}& {\delta }_{r_1}-{\nu }_1& -{\gamma }_{21}& \dots & -{\gamma }_{k1}\\ 0& -{\beta }_{v{r}_2}& -{\gamma }_{12}& {\delta }_{r_2}-{\nu }_2& \ddots & -{\gamma }_{k2}\\ ⋮& & \ddots & \ddots & \ddots & ⋮\\ 0& -{\beta }_{v{r}_k}& -{\gamma }_{1k}& -{\gamma }_{2k}& \dots & {\delta }_{r_k}-{\nu }_k\end{array}\right)$$

$$t_{21}=\frac{-det(V_{1,2})}{det(V)}=\frac{{\beta }_{hv}det(K)}{{\delta }_{h}det(V_{1,1})}$$

$${\mathcal{R}}_0=g_{11}=({\beta }_{vh})(t_{21})=\frac{{\beta }_{hv}{\beta }_{vh}det(K)}{{\delta }_{h}det(V_{1,1})}.$$

$$H=\left(\begin{array}{ccccc}0& {\beta }_{r_{1}v}& {\beta }_{r_{2}v}& \dots & {\beta }_{r_{k}v}\\ 0& 0& 0& \dots & 0\\ ⋮& & \ddots & \ddots & ⋮\\ ⋮& & \ddots & \ddots & ⋮\\ 0& 0& 0& \dots & 0\end{array}\right),W=\left(\begin{array}{ccccc}{\delta }_v& 0& 0& \dots & 0\\ -{\beta }_{v{r}_1}& {\delta }_{r_1}-{\nu }_1& -{\gamma }_{21}& \dots & -{\gamma }_{k1}\\ -{\beta }_{v{r}_2}& -{\gamma }_{12}& {\delta }_{r_2}-{\nu }_2& \dots & -{\gamma }_{k2}\\ ⋮& & \dots & \ddots & ⋮\\ -{\beta }_{v{r}_k}& -{\gamma }_{1k}& -{\gamma }_{2k}& \dots & {\delta }_{r_k}-{\nu }_k\end{array}\right).$$

$$w_{i1}^{\prime }=\frac{(-1{)}^{i+1}det(W_{1,i})}{det(W)}=\frac{(-1{)}^{i+1}det(W_{1,i})}{{\delta }_{v}det(K)},i=2,3,\dots ,k+1$$

$${\mathcal{R}}_0^r=l_{11}=\sum _{j=1}^k({\beta }_{r_{j}v})(w_{(j+1)1}^{\prime })=\sum _{j=1}^k({\beta }_{r_{j}v})(\frac{(-1{)}^{j}det(W_{1,j+1})}{{\delta }_{v}det(K)})$$

$$det(V_{1,1})={\delta }_{v}det(K)+\sum _{j=1}^k(-{\beta }_{r_{j}v})(-1{)}^{j+1}det(W_{1,j+1})={\delta }_{v}det(K)(1-{\mathcal{R}}_0^r)$$

$${\mathcal{R}}_0=\frac{{\beta }_{hv}{\beta }_{vh}det(K)}{{\delta }_{h}det(V_{1,1})}=\frac{{\beta }_{hv}{\beta }_{vh}det(K)}{{\delta }_h{\delta }_{v}det(K)(1-{\mathcal{R}}_0^r)}=\frac{{\mathcal{R}}_0^h}{1-{\mathcal{R}}_0^r}$$

$${\mathcal{R}}_0^r=\sum _{j=1}^k({\beta }_{r_{j}v})(\frac{{\beta }_{v{r}_j}{\prod }_{i\ne j}({\delta }_{r_i}-{\nu }_i)}{{\delta }_{v}det(K)})=\sum _{j=1}^k\frac{{\beta }_{r_{j}v}{\beta }_{v{r}_j}}{{\delta }_v({\delta }_{r_j}-{\nu }_j)}=\sum _{j=1}^k{\mathcal{R}}_0^{r_j}$$

Full text

1 Introduction

Some tropical diseases are amplified by one or several reservoirs. This is the case in diseases such as Chagas disease and Leishmaniasis. Indeed, Chagas disease has a domiciliary cycle, where domestic animals act as reservoirs, and a sylvatic cycle, where mammals like rodents are reservoirs [1]. Regarding Leishmaniasis, the main reservoirs of the disease in countries of South America are dogs, but other mammals could also act as reservoirs. In this paper, we are interested in diseases that have a network of reservoirs. We are also interested in representing those diseases in a simple mathematical model where we can measure the amplification effects of the reservoirs through the basic reproductive number.

In mathematical models of infectious diseases based on ordinary differential equations, the basic reproductive number of the disease is frequently obtained using the method of the Next Generation Matrix (NGM) presented in [3]. Different interpretations of the NGM can lead to different basic reproductive numbers. In Section 5 we present the construction of the NGM that it is used in this work.

As an example, we consider the system represented in Figure 1 and equations (1). This system represents the transmission of a disease between vectors $V$ and humans $H$ with transmission rates $\beta_{vh}$, $\beta_{hv}$ and mortality rates $\delta_{v}$ and $\delta_{h}$. The disease can also be transmitted among humans with transmission rate $\nu_{h}$. This model assumes that both populations are constant, so the model is determined by the equations of the infectious populations in (1).

[TABLE]

The basic reproductive number that is obtained using the NGM depends on the interpretation of which infections are considered as new. In the system presented above we could defined human and vector infections as new infections, or only human or vector infections as new. From these three interpretation we get three basic reproductive numbers (see Subsection 5.0.1). From the Theorem 2 that is proven in [3], these three numbers are greater than one (in this case the disease free equilibrium is locally asymptotically stable), or the three numbers are less than one (in this case the disease free equilibrium is unstable). In consequence, to check the stability of a possible endemic state of a system we could take an appropriate interpretation of NGM guided by the simplicity of the calculations.

In Section 2 we propose an epidemiological model of a vector-borne disease that has a network of reservoirs that infect one another. In Section 3 we show the basic reproductive number of the simplified system (omitting infections between different reservoirs) in terms of the basic reproductive number of the human cycle and the reservoirs cycles. We also present an application to Chagas disease based on data in Colombia taken from [2]. It is shown that the disease is getting extinct as long as synergistic control is made in the number of vectors and reservoirs. In Section 4 we present the discussion and conclusions of the results presented in the Section 3. In Section 5 we present the method of the NGM and the mathematical justification of results in Section 3.

2 The model

We propose a mathematical model of a vector-borne disease that has a network of $k$ reservoirs. The state variables of the system are the Human population ($H$), the vector ($V$) and the reservoirs ($R_{i},i=1,\ldots,k$). We suppose that all the populations are constant ($N_{h}$ humans, $N_{v}$ vectors and $N_{r_{i}}$ reservoirs of the species $R_{i}$, $i=1,\ldots,k$). We assume that in each reservoir species there could be self infection. Besides, the reservoirs can infect one another but there is no infection between reservoirs and human as the lines in Figure 2 shows. The parameters of the model are presented in Table 1 and the system of differential equations for the infectious populations of humans, reservoirs and vectors ($I_{h},I_{r_{i}},I_{v}$ respectively) that describes the model is given in (2).

[TABLE]

3 Results

3.1 Basic reproductive number of simplified model

For the graph in Figure 2 we define the weight of each edge as the expression next to it. For example, the weight of the edge from $H$ to $V$ is $\beta_{hv}/\delta_{h}$. The meaning of the weight of the edge from node $x$ to node $y$ is the number of infections in node $y$ that one individual of node $x$ can cause during its generation. For a cycle of the graph, we say that its weight is the geometric mean of the weights of its edges. For example, the two nodes cycle formed by $H$ and $V$ has a weight $\sqrt{\beta_{vh}/\delta_{v}\beta_{hv}/\delta_{h}}$. We denote the weight of the cycle with nodes $x$ and $y$ by $\mu_{xy}$. The following result due to Friedland [4, Theorem 2] gives us upper and lower bounds of the basic reproductive number in terms of the weights of the cycles of the graph.

Theorem 1

Let $A$ be a matrix with nonnegative entries and

[TABLE]

For $\sigma\in\Omega$, we define $\mu_{\sigma}:=(\prod_{j=i}^{k}a_{i_{j},\sigma(i_{j})})^{1/k}$. If $\mu_{*}(A):=max_{\sigma\in\Sigma}\{\mu_{\sigma}\}$ and $S(A):=(Sign(a_{i,j}))$, then

[TABLE]

and

[TABLE]

If $G$ is the NGM of an epidemiological model, then $S(G)$ determines the pairs of species where there is infection. Moreover, if $\mu_{*}$ is the heaviest cycle, we get that:

[TABLE]

In particular, this shows that a cycle with node $x$ and $y$ has basic reproductive number greater than one, i.e., $\mu_{xy}>1$, then the basic reproductive number of the whole system is also greater than one.

Constructing the NGM of the model presented in Section 2 as it is explained in Subsection 5.0.1 of Appendix, we obtain that if we consider the infection of all species as new, the matrices $F$ and $V$ that define the NGM are:

[TABLE]

In consequence, the NGM of the system is:

[TABLE]

Let us consider the system presented in the previous section with $\gamma_{ij}=0$ for all $i,j\in\{1,\ldots,k\},i\neq j$. In this case, the spectral radius of the matrix $G$ in (4) is the greatest root of the equation in (5) for $\lambda$.

[TABLE]

In general, equation (5) is not easy to solve. However, if we omit self infection in all reservoirs, i.e., $\nu_{v}=0,\nu_{h}=0,\nu_{r_{i}}=0$ for $i=1,\ldots,k$, we get that the greatest solution of (5) is given by (6).

[TABLE]

In this scenario, the set of cycles $\Omega$ would only have two nodes cycles. Moreover, If $\mu_{*}:=max_{\sigma\in\Omega}\mu_{\sigma}=max\{\mu_{vh},\mu_{vr_{1}},\ldots,\mu_{vr_{k}}\}$, the inequalities in (3) would give us the obvious bounds in (7).

[TABLE]

If we take into account the self infection in all reservoirs, the inequalities in (3) would turn into the inequalities in 8 if

[TABLE]

[TABLE]

If $\mu_{*}>1$, the disease free equilibrium would be unstable (see Theorem 2 in Appendix). If $(\sqrt{k+1}+1)\mu_{*}<1$, the disease free equilibrium would be locally asymptotically stable. Nonetheless, the inequalities in (8) does not let us determine whether $\mathcal{R}_{0}<1$ or $\mathcal{R}_{0}>1$ when $(\sqrt{k+1}+1)\mu_{*}>1$ and $\mu_{*}<1$. To solve this problem we interpret the NGM to obtain matrices $F$ and $V$ that ease the computation of the spectral radius of $FV^{-1}$.

For simplicity of the explanation, let us consider the model presented in Section 2 with only one reservoir, as the graph in Figure 3 represents.

We define the threshold values in Table 2 from their respective interpretation of the next generation matrix.

As it is shown in the Subsection 5.1, we obtain the equation (9).

[TABLE]

In consequence, $\mathcal{R}_{0}>1$ if and only if $\mathcal{R}_{0}^{h}+\mathcal{R}_{0}^{r}>1$. Using this equivalence we could determine whether $\mathcal{R}_{0}<1$ or $\mathcal{R}_{0}>1$ based on the values $\mathcal{R}_{0}^{h},\mathcal{R}_{0}^{r}$.

In Figure 4 we fix $\mathcal{R}_{0}^{h}=0.4$ and for different values of $\mathcal{R}_{0}^{r}$ we plot the stable points of the three infectious populations. In this figure we find a bifurcation in $1-\mathcal{R}_{0}^{h}=0.6$. In this example we observe how the weights of the cycles could be small but using $\mathcal{R}_{0}^{h},\mathcal{R}_{0}^{r}$ we can determine whether $\mathcal{R}_{0}<1$ or $\mathcal{R}_{0}>1$.

In the general scenario of the model presented in Section 2, we can obtain the same result in equation (9) using the values defined in Table 3.

As it is shown in the Subsection 5.1, we obtain the equation (10). If we also have that $\gamma_{ij}=0$, $i,j\in\{1,\ldots,k\}$, we obtain the equation (11).

[TABLE]

[TABLE]

In consequence, $\mathcal{R}_{0}>1$ if and only if $\mathcal{R}_{0}^{h}+\mathcal{R}_{0}^{r}>1$. Furthermore, if $\gamma_{ij}=0$, $i,j\in\{1,\ldots,k\},i\neq j,$ we get that $\mathcal{R}_{0}>1$ if and only if $\mathcal{R}_{0}^{h}+\sum_{j=1}^{k}\mathcal{R}_{0}^{r_{j}}>1$. This equivalence lets us determine whether $\mathcal{R}_{0}<1$ or $\mathcal{R}_{0}>1$ for small cycle weights, improving the informations obtained from the inequality in 8.

3.2 Application to Chagas disease

From [2], we can take some parameters for Chagas disease in Table 4. That paper considers a model with two kind of non-human host; the domiciliary hosts $R_{1}$ and the sylvatic hosts $R_{2}$. We consider the model presented in Section 2 with two reservoirs where there is no self infection ($\nu_{1}=0,\nu_{2}=0$) and there is no transmission between reservoirs ($\gamma_{12}=0,\gamma_{21}=0$). Figure 5 is the graph of this model.

If we define $\mathcal{R}_{0}^{h},\mathcal{R}_{0}^{r_{1}},\mathcal{R}_{0}^{r_{2}}$ as in Table 3, from equations (10) and (11) we have that $\mathcal{R}_{0}>1$ if and only if $\mathcal{R}_{0}^{h}+\mathcal{R}_{0}^{r_{1}}+\mathcal{R}_{0}^{r_{2}}>1$. Let us define $\mu^{\prime}_{vh}:=\frac{\beta^{\prime}_{vh}\beta^{\prime}_{hv}}{\delta_{v}\delta_{h}}$, $\mu^{\prime}_{vr_{1}}:=\frac{\beta^{\prime}_{vr_{1}}\beta^{\prime}_{r_{1}v}}{\delta_{v}\delta_{r_{1}}}$ and $\mu^{\prime}_{vr_{2}}:=\frac{\beta^{\prime}_{vr_{2}}\beta^{\prime}_{r_{2}v}}{\delta_{v}\delta_{r_{2}}}$. We have that $\mathcal{R}_{0}^{h}+\mathcal{R}_{0}^{r_{1}}+\mathcal{R}_{0}^{r_{2}}=N_{v}N_{h}\mu^{\prime}_{vh}+N_{v}N_{r_{1}}\mu^{\prime}_{vr_{1}}+N_{v}N_{r_{2}}\mu^{\prime}_{vr_{2}}$. Using the parameters in Table 4 we obtain $\mu^{\prime}_{vh}=\frac{2}{3}(\beta^{\prime}_{hv})^{2}$, $\mu^{\prime}_{vr_{1}}=\frac{2}{5}(\beta^{\prime}_{hv})^{2}$ and $\mu^{\prime}_{vr_{2}}=\frac{2}{3}(\beta^{\prime}_{hv})^{2}$. In consequence, $\mathcal{R}_{0}>1$ if and only if

[TABLE]

We must remark that $\mathcal{R}_{0}$ is always greater than $1$ if $\frac{1}{N_{v}\mu^{\prime}_{vr_{2}}}-\frac{N_{h}\mu^{\prime}_{vh}}{\mu^{\prime}_{vr_{2}}}$ is negative. This is telling that if we would want to attack the disease, we first must control the vector. In the case that $\frac{1}{N_{v}\mu^{\prime}_{vr_{2}}}-\frac{N_{h}\mu^{\prime}_{vh}}{\mu^{\prime}_{vr_{2}}}>0$, the number of reservoirs determine whether $\mathcal{R}_{0}>1$ according to the inequality in (12), as Figure 6 shows.

4 Conclusions

Based on the model of Section 2, the endemicity of the disease in one reservoir could entail the endemicity of the disease in human population. We also conclude that human endemicity of a disease in our model could not be only explained considering the dynamics of the infection within an specific system of hosts. In an specific system, we could get a small basic reproductive number that does not explain the endemicity of the disease. In our model, we observe how the basic reproductive numbers of the cycles between each reservoir and the vector could be less than one separately. However, the sum of the effects of the reservoirs can lead to endemicity of the disease in all species. In consequence, we conclude that a large enough system of hosts that contribute to spread the infection must be identified to get rid of the endemicity of the disease. As an example of control of a disease, using the data of Chagas disease we conclude that only dropping the abundancy of the reservoirs can not extinguish the disease. The abundancy of the vectors must be dropped under certain threshold for the intervention of the reservoirs to work.

5 Appendix

5.0.1 Next generation matrix

The basic reproductive number of an infectious disease $\mathcal{R}_{0}$ can be defined as the expected number of secondary cases produced in a susceptible population that are caused by an infectious individual. The NGM method lets us compute $\mathcal{R}_{0}$ in an epidemiological model where the individuals are classified in different compartments and the dynamics of the size of those compartments is described by a system of ordinary differential equations (this method is explained in [3]). The number that we get is a threshold for the local asymptotic stability of the disease free equilibrium $x_{0}$.

Let us assume that we have $N=n+m$ types of individuals and that $\bar{x}=(x_{1},\ldots,x_{n})$ represents the number of individuals in infectious compartments $I_{1},\ldots,I_{n}$ and $\bar{y}=(y_{1},\ldots,y_{m})$ represents the number of individuals in non-infectious compartments $J_{1},\ldots,J_{m}$. In the model presented in Section 2 we assume that all species populations are constant, so the system is only determined by the equations of the infectious compartments. For simplicity, we omit the equations for non-infectious compartments in this explanation. For a general exposition of the NGM, see [3].

Let us assume that the number of infectious individuals follows the system of equations in (13). In this system $f_{i}$ represents the new infections rate for the compartment $I_{i}$ and $v_{i}$ represents the rate of change of the size of this compartment due to other reasons, such as recovery, death or movement from other compartment due to causes different from new infection, like a disease stage. The functions $f_{i}$ and $v_{i}$ depend on the interpretation of which infectious individuals are regarded as new infections. Different interpretations will lead us to different versions of NGM.

[TABLE]

We define the disease free equilibrium (DFE) as an equilibrium $x_{0}=(\bar{x_{0}},\bar{y_{0}})$ of the system in (13) where $\bar{x}=(0,\ldots,0)$. We define an endemic equilibrium as an equilibrium $x_{0}=(\bar{x_{0}},\bar{y_{0}})$ where $\bar{x}\neq(0,\ldots,0)$.

Let us assume that $f_{ij}$ is the number of infections in the compartment $I_{i}$ that are caused by an individual in the compartment $I_{j}$ per unit of time in a susceptible population. In terms of the system in (13), we get $f_{ij}=\partial f_{i}(\bar{x},\bar{y})/\partial x_{j}|_{(\bar{x},\bar{y})=x_{0}}$. Let $F$ be the matrix $(f_{ij})_{n\times n}$. Let us also assume that $t_{jk}$ is the time that an individual from the compartment $I_{k}$ will be in the compartment $I_{j}$. It turns out that $V^{-1}=(t_{ij})_{n\times n}$, where $V:=(v_{ij})_{n\times n}$ and $v_{ij}=\partial v_{i}(\bar{x},\bar{y})/\partial x_{j}|_{(\bar{x},\bar{y})=x_{0}}$. If $g_{ik}$ is the number of infections in the compartment $I_{i}$ caused by an individual in the compartment in $I_{k}$ during its generation in a susceptible population, we should have $g_{ik}=f_{i1}t_{1k}+\ldots+f_{in}t_{nk}$. Each term $f_{ij}t_{jk}$ accounts for the infections caused by the individual that started in the compartment $I_{k}$ and spent a time $t_{jk}$ in the compartment $I_{j}$. Let $G$ be the matrix $(g_{ij})_{n\times n}=FV^{-1}$. We call $G$ the next generation matrix of the system in (13) with its respective interpretation contained in the functions $f_{i},v_{i},i=1,\ldots,n$. Finally, we define the basic reproductive number $\mathcal{R}_{0}$ (which depends on $G$) as the spectral radius of $G$, i.e., $\mathcal{R}_{0}=\rho(FV^{-1})$.

If the functions $f_{i}$ and $v_{i}$ satisfy the axioms $A1-A5$ presented in

[3], we obtain Theorem 2

([3, Theorem 2]).

Theorem 2

Let $x_{0}$ be the DFE of (13). Then, $\mathcal{R}_{0}<1$ implies $x_{0}$ is locally asymptotically stable and $\mathcal{R}_{0}>1$ implies that $x_{0}$ es unstable.

As an example, we consider the system of Figure 1 and equations (1) described in the Introduction. We obtain the basic reproductive numbers in Table 5.

Using Theorem 2 we have that $\mathcal{R}_{0}^{1}>1\Longleftrightarrow\mathcal{R}_{0}^{2}>1\Longleftrightarrow\mathcal{R}_{0}^{3}>1\Longleftrightarrow$ $DFE$ is unstable. In consequence, in order to verify the possible endemicity of the disease we could consider interpretations that simplify calculations.

5.1 Basic reproductive numbers of the model

For the general model defined in Section 2 we consider the numbers $\mathcal{R}_{0}$, $\mathcal{R}_{0}^{h}$, $\mathcal{R}_{0}^{r_{i}}$ and $\mathcal{R}_{0}^{r}$ defined from the systems and interpretation in Table 3. We can get the equations (14) and (15) in a similar way to the numbers presented in Table 5 in the previous Subsection.

[TABLE]

[TABLE]

As we only interpret infected humans as new infections, we get the NGM for $\mathcal{R}_{0}$ through the following matrices $F$ and $V$:

[TABLE]

Let us prove that $\rho(FV^{-1})=\frac{\mathcal{R}_{0}^{h}}{1-\mathcal{R}_{0}^{r}}$. The adjugate matrix of $V$ enables us to obtain $V^{-1}=(t_{ij})_{(k+2)\times(k+2)}$. In particular, we are interested in the entry $t_{21}$ in (16), where $V_{i,j}$ denotes the matrix that is obtained omitting the row $i$ and the column $j$ of $V$ and $K$ is the block matrix of $V$ formed by the entries $v_{ij}$ where $i=3,\ldots,k+2$ and $j=3,\ldots,k+2$.

[TABLE]

Let $G=(g_{ij})_{(k+2)\times(k+2)}$ be $FV^{-1}$. Using (16), we get that:

[TABLE]

On the other hand, when humans are not considered, we can obtain $\mathcal{R}_{0}^{r}$ from the spectral radius $\mathcal{R}_{0}^{r}=\rho(HW^{-1})$, where

[TABLE]

From the adjugate matrix of $W$, we have that if $W^{-1}=(w^{\prime}_{ij})_{(k+1)\times(k+1)}$, then

[TABLE]

Let $L=(l_{ij})_{(k+1)\times(k+1)}$ be $HW^{-1}$. Using (18), we get that:

[TABLE]

Furthermore, we have that:

[TABLE]

Replacing (20) in (17), it turns out that:

[TABLE]

If we put $\gamma_{ij}=0$, we have $det(K)=\prod_{i}(\delta_{r_{i}}-\nu_{i})$ and from equation (19) we get that:

[TABLE]