DYNAMICAL BEHAVIOR OF A ROTAVIRUS DISEASE MODEL WITH TWO STRAINS AND HOMOTYPIC PROTECTION

. A two-strain rotavirus model with vaccination and homotypic protection is proposed to study the survival of the two strains of rotavirus within the host. Corresponding to the diﬀerent eﬃcacy of monovalent vaccine against diﬀerent strains, the vaccination reproduction numbers of the two strains and the reproduction numbers of their mutual invasion are found. Based on the existence and local stability of equilibria, our results suggest that the obtained reproduction numbers determine together the dynamics of the model, and that the two-strain rotavirus dies out as both the numbers is less than unity. The coexistence of two strains, one of which is dominant, is also related to the two reproduction numbers.


1.
Introduction. Rotavirus(RV) is one of the leading causes of severe diarrhoea in children under five [4], severe cases of which can lead to fatal gastroenteritis and even death from dehydration [10]. Rotavirus infection causes 138.5 million cases of gastroenteritis and 453,000 deaths in children under 5 years of age worldwide every year, and is a serious global public health problem [14,15].
There are seven rotavirus groups. For a few groups, some different strains can be produced by recombination of their different genes [3]. This results in the broad diversity within the strains of rotavirus [8,23]. There is evidence that, for children infected by one strain, the secondary infection is more likely from another strain due to homotypic protection [18], but the corresponding incidence could decrease due to the partial immunity obtained from the infection of previous strain [6]. Immune response is also an important factor influencing infectious disease models [21]. For rotavirus infection, the immunity of children can also be obtained in two other ways. One part is the immunity inherited from the maternal antibody for infant [2,7], another acquires immunity by vaccination [1].
Vaccination programs are considered to be the most effective public health strategies for reducing the incidence rate of rotavirus infection [22,17,11,19,14,13]. Omondi, et al. [9] and Shim, et al. [12] established respectively an ordinary differential equation model and a partial differential equation model to investigate the transmission of single strain rotavirus under vaccination. So far, models including multiple strains of rotavirus infection are rare. White et al. [20] established and investigated dynamics of a series of two-strain rotavirus transmission models according to the heterogeneity of the multi-strain rotavirus for immunity. Young et al. [22] considered a more complex case that the model involves vaccination and vaccine attenuation, and examined the effects of vaccine on the qualitative behaviors of infection levels in a population. But only local dynamics of the models were obtained in [20,22] because of the complexity of the models.
Motivated by the character of rotavirus transmission mentioned above and the corresponding modelling ideas, we will propose a new dynamic model of rotavirus infection with two strains, vaccination, homotypic protection and the partial immunity. The paper is organized as follows. In next section, we present the model based on some assumptions and discuss its boundedness. In Section 3, the existence of equilibria is determined, and the corresponding regions are shown in the plane determined by the related reproduction numbers. The stability of the boundary equilibria are analyzed in Section 4. At last, the conclusion and discussion are stated.
Please use this AIMS template to prepare your tex file after the paper is accepted by an AIMS journal. Please read carefully all information including those proceeded by % sign. These are important instructions and explanations. Thank you for your cooperation.
2. Formulation and boundedness of the model. According to the introduction in the previous section, we make the following assumptions for the transmission of two strain rotavirus: 1. The infant is protected by maternal antibodies. That is, the new born has the immunity against infection of rotavirus. 2. There are two strains of rotavirus spreading in the population. Strain 2 is an evolution of strain 1, so that the children recovered from strain 2 infection could not be reinfected by strain 1. The children recovered from strain 1 infection could be reinfected by strain 2, but the incidence decreases. 3. There is vaccination in infancy. If an infant is not vaccinated, he/she will become susceptible for any strain rotavirus. 4. The vaccine here is mainly for strain 1, then the vaccinated children have complete immunity against infection of strain 1, and partial immunity for strain 2. Therefore, we divide the total population into seven subpopulations: the breast-fed infants (X), the susceptible subpopulation (S), the vaccinated subpopulation (V ), the infectious subpopulation with only strain i (I i ) (i = 1, 2), the recovered subpopulation from infection with strain i (R i ) (i = 1, 2). The transmission mechanism of the two strains of rotavirus is given by the diagram (Figure 1), where X = X(t), S = S(t), V = V (t), I i = I i (t) and R i = R i (t) (i = 1, 2) represent the numbers of individuals in the corresponding compartments at time t, respectively. A is the birth rate of the population, µ is the natural death rate of the population, τ is the rate at which maternal protection against rotavirus infection wears off, φ is the fraction of vaccination coverage for the individuals leaving the breast-fed infant subpopulation, β i (i = 1, 2) is the transmission coefficients of strain i (i = 1, 2), k is the rate at which the immunity from the vaccine weakens, ξ reflects the efficiency of the vaccine against strain 2, r i (i = 1, 2) is the rate of recovery from infection by strain i (i = 1, 2), σ reflects the homotypic protection and represents the reduction fraction of infection force of strain 2 for the individuals recovered from the infection of strain 1 relative to the susceptible individuals. Here 0 < σ, ξ, φ < 1 and other coefficients are all positive constants. Figure 1. Progression of infection from the breast-fed infants X through the susceptible S, the vaccinated V , the infected I 1 , I 2 and recovered R 1 , R 2 for the compartments of the model.
Corresponding to Figure 1, we have the following model, It is easy to know from the first equation of (1) that lim t→∞ X(t) = A τ +µ . And the variable R 2 does not appears in other equations of (1). Then the system consisting of the middle five equations of (1) has the following limit system, where M = τ A τ +µ . It is easy to know that solutions of (2) with the nonnegative initial conditions keep nonnegative. Then, from (2) it follows that It implies that lim sup From the second equation of (2), we have Hence, for an arbitrary number ε > 0, there is a T large enough such that V < V (0) + ε for t > T . Thus, for t > T , from the first equation of (2) we have It follows that Due to the arbitrariness of ε, we can get Summarizing the inferences above, we know that the set Ω = (S, V, I 1 , I 2 , R 1 ) ∈ R + 5 |S + V + I 1 + is positively invariant and attractive for system (2). Therefore, hereafter we consider the dynamic behavior of (2) on Ω.
3. Existence of equilibria. Set For the boundary equilibria of system (2), we have the following statements.
Proof. By direct calculation, we get the coordinates of E 0 and E 1 . For the equilibrium E 2 S (2) , V (2) , 0, I 2 , 0 (I 2 > 0) of system (2), its coordinates are determined by the equations From the second equation of system (6) it follows that which is less than V (0) for I 2 > 0. Substituting (7) into the first equation of (6), we have which is less than S (0) for I 2 > 0. Further, substituting equations (7) and (8) into the last equation of (6) yields, It is obvious that function f (I 2 ) is decreasing, and that lim I2→∞ f (I 2 ) = − r2+µ β2 < 0. Then, equation f (I 2 ) = 0 has a unique positive root if and only if f (0) > 0, i.e., R 20 > 1. Correspondingly, system (2) has a unique strain 2 endemic equilibrium E 2 if and only if R 20 > 1.
On the other hand, substituting (13) into the fourth equation of (11), it follows that Again, substituting it into the last equation of (11), we have Obviously, that θ 2 > θ 1 is necessary for f 2 (I 2 ) > 0. Therefore, in the following, we consider two cases: is increasing and keeps positive for I 2 > 0, and lim I2→∞ f 2 (I 2 ) = +∞. Hence, under the condition that R 10 > 1 and has a unique positive solution when R 10 > 1 and From the definition of R 10 and R 20 it follows that Then substituting them into (14), we have Then Thus, for the case that R 10 > 1 and that It is equivalent to H(θ 1 , θ 2 ) > 0, where Since which is determined by equation

KUN LU, WENDI WANG AND JIANQUAN LI
Note that the inequalities (16) and (21) can be merged into the following inequalities: For simplicity, we denote Then, from the discussions on the cases (i) and (ii), we know that, when R 10 > 1, equation 20 hold. Therefore, we have the following results with respect to the existence of positive equilibrium of system (2).
On the other hand, we define as R 10 > 1, and as R 20 > 1, where S (1) , S (2) , V (1) and R (1) 1 are defined in (5). Then, with respect to R 21 and R 12 , we have the following statements.

ROTAVIRUS DISEASE MODEL WITH TWO STRAINS
R 20 Figure 2. The existence of equilibria of system (2).
Proof. (i) When R 10 > 1, from (5) and the definition of R 10 it follows that Further, applying the expressions of R 20 and R 20 when R 10 > 1. (ii) Note that S (2) and V (2) satisfy β 2 S + (1 − ξ)β 2 V = r 2 + µ as R 20 > 1. Then applying the definition of V (0) gives 2 +k+µ < θ 2 − θ 1 , that is, According to the definition of I 2 >Î 2 is equivalent to f (Î 2 ) > 0, that is, which is the same as the inequality (18). From the proof of the existence of positive equilibrium of (2) it follows that the inequality (26) holds if and only if R 20 < R 20 . Therefore, when R 20 > 1, R 12 > 1 is equivalent to R 20 < R 20 . The proof of Proposition 1 is complete.
Thus, according to Proposition 1 Theorem 3.2 can be restated as follows. Theorem 3.3. System (2) has a unique positive equilibrium E * if R 10 > 1, R 20 > 1, R 12 > 1, and R 21 > 1, i.e., min{R 10 , R 20 , R 12 , R 21 } > 1. 4. Stability of the boundary equilibria of system (2). In this section, we first investigate the local stability of the boundary equilibria of system (2), and then discuss their global stability. The necessary and sufficient conditions for the local stability of all boundary equilibria and the global stability of disease-free equilibrium are obtained. But for the boundary equilibria E 1 and E 2 we have only the sufficient conditions for global stability.

4.1.
Local stability of the boundary equilibria of system (2). With respect to the local stability of the boundary equilibria of system (2), we have the following results.
Theorem 4.1. (i) The disease-free equilibrium E 0 of system (2) is locally asymptotically stable if and only if R 10 < 1 and R 20 < 1, unstable if either R 10 > 1 or R 20 > 1.
Proof. (i) Direct calculation shows that the eigenvalues of the Jacobian matrix of system (2) When R 10 < 1 and R 20 < 1, all the eigenvalues are negative, E 0 is locally asymptotically stable. If R 10 > 1 or R 20 > 1, the disease-free equilibrium E 0 is unstable.

ROTAVIRUS DISEASE MODEL WITH TWO STRAINS
3383 are its two eigenvalues, and the the other three eigenvalues are determined by the matrixJ We claim that all the eigenvalues ofJ(E 2 ) are with negative real parts. In fact, direct calculation shows that the characteristic equation ofJ(E 2 ) is λ 3 + c 1 λ 2 + c 2 λ + c 3 = 0, where 2 > 0. Then (2) to the last expression yields, we have 2 + k + µ > 0. So the claim holds by the Routh-Hurwitz Criterion. Hence, the local stability of E 2 is determined by λ 2 = (r 1 + µ)(R 12 − 1). Therefore, Theorem 4.1(iii) holds.
According to the results on the existence of equilibria of system (2) in Section 3, we show the local stability of its boundary equilibria in Figure 3, which corresponds to Figure 2

4.2.
Global stability of the boundary equilibria of system (2). In the following, we investigate the global stability of the boundary equilibria of system (2) by applying the theory of limit systems and constructing the appropriate Lyapunov functions. In order to make the proof of the global stabilities of equilibria E 0 , E 1 and E 2 simple, we first give the following statements about the change of the variables I 1 , R 1 and I 2 .
With respect to the global stability of the boundary equilibria, E 0 , E 1 and E 2 of system (2), we have the following results.  (2), if R 10 < 1 and R 20 < 1, the disease-free equilibrium E 0 is globally asymptotically stable on Ω.
(ii) When R 10 > 1, the strain 1 endemic equilibrium E 1 is globally asymptotically stable in Ω if one of the two sets of conditions in Lemma 4.3 holds.
Proof. (i) When R 10 < 1, by Lemma 4.2 we have that lim t→∞ I 1 (t) = lim t→∞ R 1 (t) = 0. Then, for this case, the fourth equation of (2) has the limit equation For S ≤ S (0) and V ≤ V (0) , from (33) we have So lim t→∞ I 2 (t) = 0 for R 20 < 1. Thus, the equation is the limit equation of the second equation of (2). Then lim t→∞ V = V (0) . Lastly, from lim t→∞ I 1 (t) = lim t→∞ I 2 (t) = 0 and lim t→∞ V (t) = V (0) obtained above for R 10 < 1 and R 20 < 1, the first equation of (2) has the limit system It implies that lim t→∞ S(t) = S (0) . The above inferences show that, as R 10 < 1 and R 20 < 1, equilibrium E 0 is globally asymptotically attractive. Therefore, the local stability implies that E 0 is globally asymptotically stable on Ω.
(ii) Under the condition (d1) 1 − ξ − σ ≥ 0 and R 20 < 1 or (d2) 1 − ξ − σ < 0 and R 20 ≤ min {ρ 1 , ρ 2 }, Lemma 4.3 implies that the system, consisting of the first three equations of (2), has the following limit system By the expressions of coordinates of equilibrium E 1 , (37) can be rewritten as Define a Lyapunov function L 1 by Then the derivative of L 1 with respect to t along solutions of (38) is given by The property that the arithmetic mean is greater than or equal to the geometric mean implies that dL1 dt ≤ 0, and that the equality holds if and only if S = S (1) and V = V (1) . From the first equation of (38) we know easily that singleton is the largest invariant set of system (38) (i.e. (37)) on the set satisfying dL1 dt = 0. By LaSalle Invariance Principle [5] Theorem 4.4 (ii) holds. (iii) By Lemma 4.2, R 10 < 1 implies that lim t→∞ I 1 (t) = lim t→∞ R 1 (t) = 0. Then, when R 10 < 1, the system, consisting of the first three equations of (2), has the limit system Since S (2) , V (2) and I (2) 2 satisfy the equations then system (39) can be rewritten as Define a Lyapunov function L 2 by Then the derivative of L 2 with respect to t along solutions of (40) (i.e. (39)) is (2) . By the inequality of arithmetic and geometric means (AM-GM inequality) we know that dL2 dt ≤ 0 and that dL2 dt = 0 if and only if S = S (2) and V = V (2) . Further, it is easy to see that the largest invariant set of (40) (i.e. (39)) in the set satisfying dL2 dt = 0 is the singleton Ē 2 S (2) , V (2) , I (2) 2 . Hence, according to the local stability of equilibrium E 2 of system (2), LaSalle Invariance Principle [5] implies that equilibriumĒ 2 of (40) (i.e. (39)) is globally asymptotically stable. Therefore, from the theory of limit system [16] it follows that equilibrium E 2 of system (2) is globally asymptotically stable in Ω.
According to Theorems 4.1 and 4.4, The local stability of E 0 implies its global stability. Then the analysis for the stability of E 0 is complete. When R 10 < 1 < R 20 , the local stability of E 2 also implies its global stability. For the case that R 10 ≥ 1 and R 20 > R (2) 20 , the numerical simulation shows that E 2 is also globally stable (Figure 4) although the global stability of E 2 is not proved.
Clearly, Theorem 4.4 provides only sufficient conditions for the global stability of E 1 and no condition for E * . However, the numerical simulations show that the local stability of E 1 can imply its global stability ( Figures 5 and 6), and that E * is globally stable only if it exists (Figure 7).

5.
Conclusion and discussion. For model (2) proposed here, we obtained almost complete results on the existence and local stability of equilibria of the model except for some critical cases (such as R 10 = 1 or R 20 = 1). With respect to the corresponding conditions, R 10 and R 20 are two important quantities. Note that, in the presence of vaccination which is fully effective for the strain 1 and partially effective for the strain 2, S (0) and S (0) + (1 − ξ)V (0) are the numbers of susceptible individuals for the two strains without infection, respectively. And 1 ri+µ (i = 1, 2) is the average infectious period of the strain i. Therefore, is the vaccination reproduction number of the strain 1 rotavirus for model (2), and is the vaccination reproduction number of the strain 2 rotavirus for model (2).  Figure 6. The trajectories of I 1 and I 2 for the case that R 10 > 1 and 1 < R 20 < R 20 . Here, β 1 = 2 and β 2 = 0.2.
Correspondingly, R 10 = 44.28, R 20 = 4.54, R (1) 20 = 5.26, E 1 is globally stable. In addition, we also give two quantities R 12 and R 21 , which are defined in (24) and (25), respectively. As shown in Theorems 3.3 and 4.4, they determine the existence of the positive equilibrium E * and the stability of equilibria E 1 and E 2 .
The condition that R 21 is defined is that R 10 > 1. For this case, E 1 exists. According to assumption that the children recovered from strain 1 infection could be reinfected by strain 2, then when population stay in equilibrium E 1 , the individuals susceptible for strain 2 consist of S, V and R 1 . Therefore, by the biological meanings of the parameters ξ and σ we know easily that R 21 is the invasion reproduction number of strain 2 rotavirus under the case that the infection of strain 1 is at the steady state.
The condition that R 12 is defined is that R 20 > 1, which assures that E 2 exists. According to assumption that the vaccinated children have complete immunity against infection of strain 1 and that the children recovered from strain 2 infection could not be reinfected by strain 1, then when population stay in equilibrium E 2 , the individuals susceptible for strain 2 include only S. Therefore, R 12 is the invasion reproduction number of strain 1 rotavirus under the case that the infection of strain 2 is at the steady state.
Thus, Theorem 3.3 suggests that the coexistence of the two strains implies the success of their mutual invasion.
With respect to the global stability of the boundary equilibria of the model, the almost necessary and sufficient conditions are found for the infection-free equilibrium E 0 by applying the theory of limit systems and the comparison theorem, but only the sufficient conditions are obtained for the single strain equilibrium by applying the theory of limit systems and constructing the appropriate Lyapunov functions. Many numerical simulation showed that the local stability of the strain i (i = 1, 2) equilibrium implies its global stability.
For the positive equilibrium of the model (coexisting equilibrium of the two strains), E * , only the existence is determined, the condition on its local and global stability are not obtained. However, our numerical simulation shows that E * should be globally stable if it exists. For these results shown by numerical simulation but not obtained by theoretical analysis, we will continue to explore the analytical methods in subsequent studies.