Multiple infection leads to backward bifurcation for a schistosomiasis model

1.   Introduction

Schistosomiasis is one of neglected tropical diseases. On February 25th, 2016, WHO reported that at least 258 million people required preventive treatment for schistosomiasis in 2014 ^[16]. More than 61.6 million people were reported to have been treated for schistosomiasis in 2014. People can be infected by being exposed to infested water during routine agricultural, domestic, occupational and recreational activities. Especially in the flood period, many people may be infected multiply by infested water ^[9,18]. Hence, schistosoma infections remain a serious public health problem worldwide. Identification of critical factors in the spread of schistosomiasis can be helpful for the control and prevention of the disease.

Although drug therapies for treating schistosomiasis are available, multiple infection is still a serious problem, as it increases the difficulty for schistosomiasis control. Xu et al. reported that more than 90 percent of the schistosomiasis patients were repeatedly exposed to infected water ^[17]. The study of Zhao et al. ^[19] involves a total of 679 people, of which 175 people are uninfected, 226 people were infected once, 106 people were infected 4 times and 172 people were infected 6 times. Different number of exposures may lead to various degrees of infection. It is reported in ^[1] that the infectivity of an infected human host is dependent on the infection history. Hence, multiple infection may cause the patient's condition is more serious.

Studies about treatment-reinfection were done in a region of Philippines where S. japonicum is endemic. In these studies, the impact of reinfection with S. japonicum after treatment with praziquantel on the mean hemoglobin level is evaluated. Results show that rapidly reinfected individuals did not have the positive treatment effect on hemoglobin seen in non-reinfected individuals ^[11].

From the perspective of epidemiology, multiple infection may greatly influence the transmission and control of schistosomiasis. Hence, it is important to increase our understanding of multiple infection in schistosomiasis transmission dynamics, which may help identify more effective measures in containing or eliminating its transmission. Currently, many mathematical models have been developed for investigating schistosomiasis dynamics ^{[5,7,8,13,14,20]} These models have provided useful information for understanding the mechanics of schistosomiasis transmission.

However, these models use the same threshold quantity to determine whether or not schistosomiasis is endemic. The threshold quantity is the basic or effective reproduction number. If the reproduction number is greater than 1, schistosomiasis is persistent and dies out if it is below 1. In this paper, multiple infection of patients was incorporated in the mathematical model. It is found that a backward bifurcation may occur if the infectivity of people with multiple infection is sufficiently high. This suggests the importance of studying the effects of multiple infection on the spread and control of schistosomiasis. Moreover, such results can help to obtain thresholds for the control of schistosomiasis in the case when the threshold value of the reproduction number is no longer 1.

The paper is organized as follows. In Section 2 we establish a schistosomiasis model with multiple infection of human. The dynamics for the model are studied in Section 3. This section includes the threshold condition, the existence of equilibria, backward bifurcation analysis and the global stability of the disease free equilibrium. Some numerical simulations are presented in Section 4.

2.   Modeling

Based on years of experience in schistosomiasis prevention and treatment, one of the typical features of schistosomiasis is multiple infection of a human host by parasites, which may dramatically affect the host's infectivity. The problem is whether multiple infection may cause greater damage than single infection or not ^[11,17]. Therefore, it is necessary to study the effects of multiple infection on the spread of schistosomiasis in theory. In this model, one classified human three classes, that is susceptible $S_{H}$, single infected $I_{H}$ and multiple infected $M_{H}$ according to the epidemiology of schistosomiasis ^[17], respectively. Similarly based on the epidemiology of snail, snails are only separated two classes: susceptible snail $S_{S}$ and infectious snail $I_{S}$.

The annual 4-10 month is schistosomiasis epidemic period, and in this period, people contact infected water with high frequency. For example, flooding occurred around July ^[9,18]. Due to the good effect of drug control, few people died of schistosomiasis ^[12]. In additional, this epidemic period is short. Hence, mortality due to illness for human and snail is ignored in this model. The death rate of three kinds of people is the same $\mu_{H}$. For human, both single infected and multiple infected can recovery with different rate due to drug treatment ^[12]. Reference [11 showed that the treatment effect of multiple infection is not good with single infection. According to the work experience of the schistosomiasis agency for many years, with the same treatment rate, people with multiple infections $M_{H}$ will enter the lightly infections group $I_{H}$ after treatment, while those with single infection $I_{H}$ will enter the susceptible group $S_{H}$. Here, $\gamma_{H}$ represents the recovery rate. The transmission rate from single infected snail to human is $\beta_{SH}$. However, once a person is infected, people will be uncomfortable ^[17]. People may choose to take a break, which will reduce the chance of exposure to infected water. Hence, we denote $\omega$ to describe the relative transmission from snail to multiple infected human with $0 < \omega\leq 1$ ^[12]. In additional, the number of eggs per gram (EPG) in the body of multiple infected people is more than in single infected people. Thus multiple infected people has stronger ability to infect snail than single infected people. Here using $\epsilon$ to describe the relative transmission from multiple infected human to snail with $\epsilon\geq 1$ ^[12].

Then the model is as the following

All the parameters are listed in Table 1. Our aim was to study the effect of multiple infection on the dynamic behavior of schistosomiasis. In detail, it is to study the effect of the two parameters $\omega, \varepsilon$ on the dynamics of the schistosomiasis transmission model.

3.   Dynamics of the model

It is easy to see that for (2.1) all trajectories in the positive cone enter or stay inside the region

That means that $\Omega$ is a positively invariant set of (2.1).

The disease free equilibrium (DFE) $E_{0} = (\frac{\Lambda_{H}}{\mu_{H}}, 0, 0, \frac{\Lambda_{S}}{\mu_{S}}, 0)$ always exists. According to the concept of next generation matrix ^[3] and the technique developed in ^[4] for calculating the basic reproduction number, if we follow the notations in ^[4] and let

Then the basic reproduction number for model (2.1) is

Considering the convenient calculation later, we still denote

Other equilibria are obtained from (2.1) by equating the right-hand side of (2.1) to zero and solving the resulting algebraic equations. For any equilibrium $E^{*} = (S_{H}^{*}, I_{H}^{*}, M_{H}^{*}, S_{S}^{*}, I_{S}^{*})$, its coordinates satisfy the following relations

and $I_{S}^{*}$ is the root of the following equation

Here

If $R_{0}>1$, $C < 0$ and then the equation (3.1) must have a unique root $I_{S}^{*}$. This means that there exists a unique endemic equilibrium $E^{*}$ if $R_{0}>1$. Note that

Moreover, from $R_{0}\leq1$ one can get $\varepsilon^{*}\geq1+\frac{\Lambda_{H}\beta_{HS}}{\mu_{S}\mu_{H}}+\frac{\omega}{\gamma_{H}+\mu_{H}}$, that is $\varepsilon^{*}$ must be bigger than $1$ when $R_{0}\leq1$.

If $R_{0} = 1$, $C = 0$ and then the equation (3.1) have positive root under the condition that $B < 0$. Hence, the equation (3.1) has a unique positive root $I_{S}^{*} = -\frac{A}{B}$, which corresponds to the unique endemic equilibrium $E^{*}$ if and only if $\varepsilon > \varepsilon^{*}$. Otherwise, the system (2.1) has no endemic equilibrium when $\varepsilon \leq \varepsilon^{*}$.

Now we consider the case when $R_{0} < 1$. In this case, the system (2.1) may have endemic equilibria only when $\varepsilon > \varepsilon^{*}$. Let the discriminant of (3.1) be $\Delta$, then we obtain

It is easy to see that $\Delta$ is a strictly monotone increasing function of $R_{0}$ in the right half plane.

Under $\varepsilon > \varepsilon^{*}$, $\Delta = 0$ is equivalent to $B = -2\sqrt{AC}$, that is

Let $\sqrt{1-R_{0}} = p \ \ (0 < p < 1)$, we can get an equation $g(p) \triangleq a p^{2}+b p+c = 0$ about $p$. Here

From $R_{0} < 1$ and $\varepsilon > \varepsilon^{*}$, we have $c < 0$. In additional, $g(1)>0$. Hence, there exists a unique solution $p^{*} = \frac{-b+\sqrt{b^{2}-4ac}}{2a} < 1$ for the equation $g(p) = 0$. That is to say that there exists a unique $R_{0}^{*} = 1-{p^{*}}^{2} < 1$ for the equation $\Delta = 0$ in terms of $R_{0}$ under the condition that $\varepsilon > \varepsilon^{*}$. Then we can obtain the following equivalence relations:

This means that the equation (3.1) may have $0$, $1$ and $2$ roots. In summary, the result about the number of endemic equilibria is in Theorem 3.1.

Theorem 3.1. For the system (2.1), with $R_{0}^{*}$ and $\varepsilon^{*}$ defined as above, we have

(1) When $R_{0}>1$, there is a unique endemic equilibrium $E^{*}$.

(2) When $R_{0} = 1$ and $\varepsilon > \varepsilon^{*}$, there is a unique endemic equilibrium $E^{*}$.

(3) When $R_{0}\leq1$ and $\varepsilon \leq \varepsilon^{*}$, there is no endemic equilibrium.

(4) When $R_{0}^{*} < R_{0} < 1$ and $\varepsilon >\varepsilon^{*}$, there are two endemic equilibrium $E^{*}$ and $E_{*}$.

(5) When $R_{0} = R_{0}^{*}$ and $\varepsilon > \varepsilon^{*}$, $E^{*}$ and $E_{*}$ coalesce at a unique endemic equilibrium of multiplicity 2.

(6) When $R_{0} < R_{0}^{*}$ and $\varepsilon > \varepsilon^{*}$, there is no endemic equilibrium.

Here, $E^{*} = (S_{H}^{*}, I_{H}^{*}, M_{H}^{*}, S_{S}^{*}, I_{S}^{*})$ and $E_{*} = (S_{H*}, L_{H*}, H_{H*}, S_{S*}, I_{S*})$ are the corresponding equilibria, and

Remark 3.2. When $\varepsilon >\varepsilon^{*}$, this theorem shows that there exists $R_{0}^{*}$ ($0 < R_{0}^{*} < 1$) such that the system (2.1) has a unique endemic equilibrium for $R_{0} = R_{0}^{*}$, then the system has two endemic equilibrium for $R_{0}^{*} < R_{0} < 1$ and a unique endemic equilibrium for $R_{0} = 1$. This situation corresponds to a backward bifurcation which occurs at $R_{0} = 1$ (Figure 1). When $\varepsilon \leq \varepsilon^{*}$, at $R_{0} = 1$ a forward bifurcation occurs at $E_{0}$. In the prevention and treatment of schistosomiasis, special consideration should be given to the patient's multiple infection and the treatment of multiple infectious patients.

Remark 3.3. Here $\varepsilon^{*}$ need to be greater than $1+\frac{\Lambda_{H}\beta_{HS}}{\mu_{S}\mu_{H}}+\frac{\omega}{\gamma_{H}+\mu_{H}}$ at least, such that there is the possible backward bifurcation when $R_{0} < 1$. Multiple infection caused the human body to be infected by more schistosoma than single infection. The number of eggs per gram (EPG) in fecal is proportional to the number of schistosoma in the human body ^[2]. This means multiple infected patient's EPG is more than single infected patient. From the data in ^[2], 0 EPG represented uninfected or susceptible, 24-96 EPG indicated lightly infected, 120-792 EPG represented moderate infection and bigger than 816 indicated severely infected. These data indicate that EPG of patients with multiple infections have at least 5 times more than single infection patients. Hence, here $\varepsilon^{*}\geq 1+\frac{\Lambda_{H}\beta_{HS}}{\mu_{S}\mu_{H}}+\frac{\omega}{\gamma_{H}+\mu_{H}}$ is realistic.

Remark 3.4. Note that the backward bifurcation will not happen if $\omega = 0$. That means if there is not multiple infection, there will be no backward bifurcation. To control schistosomiasis, the basic reproduction number $R_{0}$ is only threshold.

Now we study the stability of equilibria. It is easy to calculate that the characteristic equation about $E_{0}$ is given by

Then the characteristic roots are $-\mu_{H}$, $-\mu_{S}$ and roots of the following equation:

It is easy to see that the real parts of all eigenvalues of $E_{0}$ are negative if and only if $R_{0} < 1$. Hence, the disease free equilibrium $E_{0}$ of the system (2.1) is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0}>1$.

In order to study the global stability of the disease free equilibrium, we use Metzler matrix theory and the technique of ^[10]. Consider systems of the following form:

where $x_{1}\in R^{n_{1}}_{+}$, $x_{2}\in R^{n_{2}}_{+}$, and $f$ and $g$ are $C^{1}$. We denote by $x = (x_{1}, x_{2})$ the state of the system and $(x^{*}_{1}, 0)$ is a DFE on the positively invariant set $\Omega \subset R^{n_{1}+n_{2}}_{+}$. Now rewritten (3.2) as:

For system (3.3), we make the following assumptions:

($\textbf{H}_{1}$): The system is defined on a positively invariant set $\Omega$ of the nonnegative orthant. The system is dissipative on $\Omega$.

($\textbf{H}_{2}$): The sub-system $\dot{x}_{1} = A_{1}(x)\cdot(x_{1}-x^{*}_{1})$ is globally asymptotically stable at the equilibrium $x^{*}_{1}$ on the canonical projection of $\Omega$ on $R^{n_{1}}_{+}$.

($\textbf{H}_{3}$): The matrix $A_{2}(x)$ is Metzler and irreducible for any given $x\in\Omega$.

($\textbf{H}_{4}$): There exists an maximum matrix $\bar{A}_{2}$, then for any $\bar{x}\in\Omega$ such that $\bar{A}_{2} = A_{2}(\bar{x})$, $\bar{x}\in R^{n_{1}}_{+}\times\{0\}$.

($\textbf{H}_{5}$): $\alpha(\bar{A}_{2})\leq 0$, i.e., the greatest real part of eigenvalues of $\bar{A}_{2}$ is nonnegative.

Now we state two lemmas due to ^[10].

Lemma 3.5. If the above hypothesis $\textbf{H}_{1}$-$\textbf{H}_{5}$ are satisfied, then the DFE is globally asymptotically stable in $\Omega$.

Lemma 3.6. If the same notations and hypothesis in Lemma 3.5 hold and if furthermore we have $\bar{A}_{2} = A_{2}(x^{*}_{1}, 0)$, the DFE is globally asymptotically stable if and only if $\alpha(\bar{A}_{2})\leq 0$.

Next we discuss the global stability of the disease free equilibrium $E_{0}$ using the above two Lemmas. From the system (2.1), we know

is a compact positively invariant absorbing set contained in the nonnegative orthant. Thus the system (2.1) is dissipative on $\Omega$ because the trajectories of (2.1) are forward bounded. Now we shall study the system (2.1) on $\Omega$.

We set for system (2.1) $x_{1} = (S_{H}, S_{S})$, $x_{2} = (I_{H}, H_{H}, I_{S})$ and $x_{1}^{*} = (\frac{\Lambda_{H}}{\mu_{H}}, \frac{\Lambda_{S}}{\mu_{S}})$. As in ^[10], we express the sub-system as $\dot{x}_{1} = A_{1}(x_{1}, 0)(x_{1}-x_{1}^{*})$ and

This is a linear system and its unique equilibrium ($\frac{\Lambda_{H}}{\mu_{H}}, \frac{\Lambda_{S}}{\mu_{S}}$) (corresponding to the DFE of (2.1)) is globally asymptotically stable, hence the assumption $\textbf{H}_{1}$ and $\textbf{H}_{2}$ are satisfied.

The matrix $A_{2}(x)$ is given by

As required by hypothesis $\textbf{H}_{3}$, for any $x \in \Omega$ the matrix $A_{2}(x)$ is irreducible.

Now let us check $\textbf{H}_{4}$. There is a maximum which is uniquely realized in $\Omega$ if $S_{H} = \frac{\Lambda_{H}}{\mu_{H}}, S_{S} = \frac{\Lambda_{S}}{\mu_{S}}$ which corresponds to the DFE. This maximum matrix is then $J_{2}$, the sub-block of the Jacobian matrix at the DFE, corresponding to the matrix $A_{2}(x)$. The matrix $J_{2}$ is given by

Therefore, we are in the situation of Lemma 3.6 where the maximum is attained at the DFE.

The hypothesis $\textbf{H}_{5}$ requires that $\alpha (J_{2})\leq 0$. Rewriting $J_{2}$ into a sub-block matrix as $J_{2}$ $=$ $\left(\begin{array}{cc} A_{1 \times 1} & B_{1 \times 2} \\ C_{2 \times 1} & D_{2 \times 2} \\ \end{array} \right).$ Since $A$ is already a Metzler stable matrix, the condition $\alpha (J_{2})\leq 0$ is equivalent to $\alpha (D-CA^{-1}B)\leq 0$ ^[10] and this last condition is equivalent to the $F_{0}$ condition

We have computed $F_{0}$ and one can see that the hypothesis $\textbf{H}_{1}$, $\textbf{H}_{2}$, $\textbf{H}_{3}$, $\textbf{H}_{4}$ and $\textbf{H}_{5}$ are all satisfied. Then by Lemma 3.6, we have proved the following Theorem 3.7.

Theorem 3.7. The disease free equilibrium $E_{0}$ of the system (2.1) is globally asymptotically stable when $R_{0} < 1$.

In additional, using the Routh-Hurwitz criterion the stability of $E^{*}$ can be obtained through complex computation.

Theorem 3.8. The unique endemic equilibrium $E^{*}$ of the system (2.1) is locally asymptotically stable when $R_{0}>1$.

4.   Numerical simulation

Based on some parameter values in the Table 1, numerical simulations are conducted in this section. Here $R_{0} = 0.98$, $R_{0}^{*} = 0.95$ and $\varepsilon^{*} = 6.82$. First we can find two endemic equilibrium when $R_{0}^{*} < R_{0} < 1$ and $\varepsilon >\varepsilon^{*}$ and backward bifurcation occurs (Figure 2).

The following figures show that the unique endemic equilibrium is stable if $R_{0}>1$ (see Figure 3, 4).

Not only that, the number of single and multiple infectious patients changes as the change of $\omega$. If the value of $\omega$ is larger, the number of patients with multiple infection is greater. Although the number of patients with single infection is lower at the beginning, eventually becomes great. Moreover, very little change in $\omega$ can lead to big change in the number of patients. This means multiple infection will lead to more and more new patients with single and multiple infections.

From Figure 4 we can see the number of single and multiple infected humans also increases when the value of $\varepsilon$ is greatly increased. If the amplitude of the value of $\varepsilon$ increase is small, the impact on the number of patients is not very large. The main factors that to affect $\varepsilon$ are feces of the patients who have repeatedly infected. Therefore, to do a good job of health facilities, especially for patients with multiple infections, is good to block the spread of schistosomiasis.

5.   Discussion

Multiple infection is a frequent phenomenon in the real life, especially in the process of the transmission of schistosomiasis which is related to water, because people can not live without water. In this paper, multiple infection of patients was considered into mathematical model. It is found that a backward bifurcation occurs if transmission rate from multiple infectious human to snail is big. In additional, if there is not multiple infection, there will be no backward bifurcation. Hence, special consideration should be given to the patient's multiple infection and the treatment of multiple infectious patients in the prevention and treatment of schistosomiasis.

Based on the simulations, we found the unique endemic equilibrium is stable when the basic reproduction number is greater than 1. Furthermore, multiple infection will lead to more and more new patients with single and multiple infections. Therefore, we should vigorously publicize the risk of multiple infections, to remind people not to contact with water, especially in the flood period. In additional, If the transformation from multiple infection to snail is bigger than $\varepsilon^{*}$, to obtain stability of the disease free equilibrium the basic reproduction number should be smaller than $R_{0}^{*}$. This means the basic reproduction number is not the unique threshold to control schistosomiasis. To eliminate schistosomiasis, we should do more work to reduce the basic reproduction number to be smaller than $R_{0}^{*}$. For example, the snail density should be controlled to be less than $\frac{R_{0}^{*}\mu_{H}\mu_{S}(\gamma_{H}+\mu_{H})}{\Lambda_{H}\beta_{SH}\beta_{HS}}$.

Eggs infected by multiple infected patients can cause serious infection to snail. Numerical simulations show that this serious infection can lead to more new patients, which may lead to schistosomiasis outbreaks. In addition, multiple infections also cause great harm to the human body. Hence, prevention of multiple infection is particularly important.

Acknowledgement

The authors greatly appreciate the anonymous referees for their valuable comments which help to improve the paper. This research is supported by National Natural Science Foundation of China (11401002), Natural Science Fund for Colleges and Universities in Anhui Province(KJ2018A0029) and Quality Engineering Project of Anhui University (ZLTS2016065).

Conflict of interest

The authors declare there is no conflict of interest.