EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 13, No. 4, 2020, 840-851 ISSN 1307-5543 – www.ejpam.com Published by New York Business Global A Delayed Vaccination Model for Rotavirus Infection Florence A. Adongo1, Lawrence O. Onyango2,∗, Job Bonyo3, G.O. Lawi4, Ogada A. Elisha5 2,5 Mathematics Department, Faculty of Science, Egerton University, Nakuru, Kenya 1,3 Mathematics Department, Faculty of Science, Maseno University, Kisumu, Kenya 4 Mathematics Department, Faculty of Science, Masinde Muliro University of Science and Technology, Kakamega, Kenya Abstract. In this work, a mathematical model for rotavirus infection incorporating delay differ- ential equations has been formulated. Stability analysis of the model has been performed. The result shows that the Disease Free Equilibrium is globally asymptotically stable and the Endemic Equilibrium undergoes a Hopf bifurcation. Numerical analysis has been performed to validate the analysis. 2020 Mathematics Subject Classifications: 34C60, 34xx, 37Nxx Key Words and Phrases: Rotavirus, Hopf Bifurcation, Global asymptotic stability, Endemic Equilibrium 1. Introduction This work is an extension of the work done by Onyango et.al [8]. In their work, they assumed that the effects of vaccines are immediate. This has not been the case since there must be time lapse. It is therefore necessary to investigate how this time lapse impact on the effectiveness of the vaccine. For the literature review, proof of existence of both disease free and endemic equilibrium, the establishment of basic reproduction number of the model, see [8]. This work is organized as follows. The model is formulated in section 2, in section 3, the model is analyzed. In section 4, the numerical simulation and discussion is performed. Conclusion and discussion is done in section 5. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v13i4.3822 Email addresses: florenceeadongo@gmail.com (F.A. Adongoi), lawionyi@gmail.com (O.L. Onyango), jobbonyo@maseno.ac.ke (J. Bonyo), achiengelisha@gmail.com (E.A. Ogada), glawi@mmust.ac.ke (G.Lawi) https://www.ejpam.com 840 c© 2020 EJPAM All rights reserved. L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 841 2. Model description and formulation. We achieve the objective of this study, by formulating a mathematical model based on a system of delay differential equation for rota-virus incorporating time delay in the effects of vaccination. The total human population size understudy, N(t), at any time is subdivided into classes such as susceptible S, infectious with rotavirus I, vaccinated V and recovered R. The total population N(t) = S(t) + V(t) + I(t) + R(t). Assuming that the mass action incidence transmission is defined by βSI, where β is the effective contact rate for disease transmission and the initial conditions are such that the parameter S, V, I, R remain non-negative for all time t ≥ 0. Stability analysis of the model is done to determine the conditions for the spread of disease in a given population. Rotavirus spreads by contact with infected faeces and might also be transmitted through faecally-contaminated; food, water and respiratory droplets [1]. Since the incubation period is very short [8], we assume that the probability of survival till the infectious state for the individuals exposed to rotavirus is unity and therefore exclude the exposure stage. The individuals infected with rotavirus include both symptomatic and asymptomatic cases because they are capable of infecting others [9]. The recovery class comprises of those who have been removed from the scene of infection by such means as infection acquired immunity and death. It is possible children can develop some level of immunity to rotavirus from maternal antibody due to breastfeeding but this immunity does not last for long hence we consider the effect of vaccination at birth and vaccination of susceptibles [8]. The human population is not assumed to be constant, since birth, immigration, emigration and death occur. Assumed a constant recruitment ρ out of which (1 − ρ)Λ is into susceptible class and ρΛ is into the vaccinated class. Susceptibles are vaccinated at the rate γ and the vaccine efficacy which has been shown to wane is assumed to take place at the of ω [6]. The parameter 0 ≤ 1−ε < 1 models the decrease in the risk of infection as a result of vaccination. Disease mortality takes place at the rate δ and recovery from infection takes place at the rate κ, it therefore natural that after a single natural infection immunity is developed, subsequent infections are less severe [8]. The population decreases due to natural deaths at a rate µ. Most vaccines take time in the body to become effective, this is because immunity has to be developed to protect the body against infection. The time the vaccine is administered and the time it becomes effective is defined as time delay denoted by τ . These parameter values are summarized in Table 1 below. The flow chart for the proposed model is shown in Fig.1 Finally, from the above definition and assumptions, the proposed mathematical model and the parameters are described below: dS dt = (1− ρ)Λ− βSI − γS + ωV − µS, dV dt = ρΛ + γS − (1− ε)βV (t− τ)I − (ω + µ)V, dI dt = βSI + (1− ε)βV (t− τ)I − (δ + κ+ µ)I, L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 842 Table 1: Parameter values Parameter Symbol Recruitment rate into susceptible (1− ρ)Λ Recruitment rate into vaccination ρΛ Vaccination rate of susceptible γ Vaccine efficacy waning rate ω Expected decrease in the risk of infection ε Rate of flow into the removed class κ Transmission rate β Natural death rate of human µ Rotavirus induced deaths δ Vaccinated individual ρ Delay time in vaccination τ dR dt = κI − µR (1) We set the initial conditions for system (1) as S(θ) = φ1(θ), V (θ) = φ2(θ), I(θ) = φ3(θ), R(θ) = φ4(θ), where φi(θ) ≥ 0, θ ∈ [−τ, 0], φ(0) > 0, for i = 1, 2, 3 · · · such that φ = (φ1, φ2, φ3, φ4) are defined in the Banach space of continuous functions mapping in the interval [−τ, 0]⇒ R4. 3. Model analysis Since the model has been analysed by Onyango et.al [8] when delay has not been incorporated, it is in order that we analyse only the effect of delay on both the disease free equilibrium and local stability of the endemic equilibrium. 3.1. The Global Stability of Disease Free Equilibrium To prove the global stability, we use Rv as derived and explained by Onyango et.al, see [8]. The global stability of equilibrium is globally asymptotically stable if Rv ≤ 1. We use the technique by Castillo Chavez [4]. We write system (1) in the form; dX dt =H(X,Z) dZ dt =G(X,Z), G(X, 0) = 0 Where X ∈ R2 denotes uninfected compartments (S,V) and Z ∈ R1 denotes infected compartment (I). The disease free equilibrium is now denoted as E0 1 = (X0, 0), X0 = ( ωΛ+(1−ρ)µΛ µ(ω+γ+µ) , (γ+µρ)Λ µ(µ+ω+γ) ) The technique stipulates that the following conditions H1 and H2 must be met to guarantee global asymptotically stability: H1: For dX dt =H(X, 0), X0 is globally asymptotically stable. H2: G(X,Z) = PZ − Ĝ(X,Z), Ĝ(X,Z) ≥ 0 for (X,Z) ∈ Ω, L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 843 Figure 1: Model Flow Chart where X0 is the disease free equilibrium, P = DzG(X, 0) is an M-matrix (the off- diagonal element of P are non-negative) and Ω is the region where the model (1) is bio- logically feasible. Theorem 1. The disease free equilibrium E0 of system (1) is globally stable if Rv < 1 and unstable whenever Rv > 1, provided that the conditions H1 and H2 above are satisfied. Proof. From system (1), we can clearly see that H(X, 0) = ( (1− ρ)Λ− (γ + µ)S + ωV ρΛ + γS − (ω + µ)V ) and G(X,Z) = PZ − Ĝ(X,Z) Differentiating the right hand side of equation 3 of system (1) with respect to I, we obtain P = βS + (1− ε)βV (t− τ)− (δ + κ+ µ). Therefore PZ = βSI + (1− ε)βV (t− τ)I − (δ + κ+ µ)I and ĜZ = PZ −GZ = βSI + (1− ε)βV (t− τ)I − (δ + κ+ µ)I − (βSI + (1− ε)βV (t− τ)I − (δ + κ+ µ)I) = 0 Since conditions H1 and H2 are satisfied, the disease free equilibrium is therefore globally asymptotically stable when Rv < 1 and unstable whenever Rv > 1. L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 844 3.2. The Local Stability of Endemic Equilibrium (E.E.) The Jacobian matrix at the endemic equilibrium E∗ can be expressed as J = a1 a2 a3 a4 a5e −λτ + a6 a7 a8 a5e −λτ a9  where, a1 = −(βI∗ + γ + µ) a2 = ω a3 = −βS∗ a4 = γ a5 = −(1− ε)βI∗ a6 = −(ω + µ) a7 = (1− ε)βV ∗(t− τ) a8 = βI∗ a9 = βS∗ + (1− ε)V ∗(t− τ)− (δ + κ+ µ) To find the characteristic equation of the linearized system (1) at steady states, we compute the eigenvalues of the following matrix∣∣∣∣∣∣ λ− a1 a2 a3 a4 λ− [a5e −λτ + a6] a7 a8 a5e −λτ λ− a9 ∣∣∣∣∣∣ = 0 this gives λ3 +M2λ 2 +M1λ+M0 + (N2λ 2 +N1λ+N0)e−λτ = 0, (2) where M2 = −(a1 + a6 + a9) M1 = (a6a9 + a1a6 + a1a9 − a2a4 − a3a8) M0 = (−a1a6a9 − a2a4a9 − a2a7a8 + a3a6a8) N2 = −a5 N1 = (a5a9 − a5a7 + a1a5) N0 = (a1a5a7 − a1a5a9 + a3a4a5 + a3a5a8). Multiplying both sides of equation (2) by eλτ , we obtain (λ3 +M2λ 2 +M1λ+M0)eλτ +N2λ 2 +N1λ+N0 = 0 (3) When τ = 0, equation (3) can be expressed as λ3 +M02λ 2 +M01λ+M00 = 0, (4) L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 845 where M02 =(M2 +N2) M01 =(M1 +N1) M00 =(M0 +N0). Based on Routh Hurwitz theorem [10], it can be concluded that all the roots of equation (4) are in the open left half plane if and only if the following condition holds: C1 : M02 > 0,M00 > 0 and M02M01 > M00. If τ = 0 and C1 holds then equation (2) is locally asymptotically stable. For τ > 0, we let λ = iw, for w > 0 be a root of equation (3). Substituting λ = iw into (2) we obtain (−iw3 −M2w 2 + iM1w +M0)((cos(wτ) + isin(wτ))−N2w 2 + iN1w +N0 = 0 (5) On separating the real and imaginary parts of equation (5) gives p1(w)cos(wτ)− p2(w)sin(wτ) = p3(w) P4(w)sin(wτ) + p5(w)cos(wτ) = p6(w), (6) where p1(w) = −M2w 2 +M0 P2(w) = M1w − w3 p3(w) = N2w 2 −N0 p4(w) = M0 −M2w 2 p5(w) = M1w − w3 p6(w) = −N1w Solving equation (6), we obtain cos(wτ) = p01(w) p00(w) sin(wτ) = p02(w) p00(w) , (7) where p00 = M2w 4 − (2M0M2 −M2 1 )w2 − w6 +M0 + 2M1w 4 p01 = (M0N2 +N0M2 −M1N1)w2 − (M2N2 −N1)w4 −N0M0 P02 = (M1N2 −N0)w3 −N2w 5 −N0M1w Squaring and adding the two equations in equation (7), we get p2 01(w) + p2 02(w)− p2 00(w) = 0 (8) L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 846 Suppose that C2: equation (8) has at least one positive root, w0, then equation (4) will definitely have pure imaginary roots ±iw0. For w0 we obtain the critical value of time delay as shown below τ0 = 1 w0 arccos { p01(w0) p00(w0) } (9) Differentiating equation (3) implicitly with respect to τ we obtain ( 3λ2 + 2M2λ+M1 ) eλτ dλ dτ + ( λ+ τ dλ dτ )( λ3 +M2λ 2 +M1λ+M0 ) eλτ+(2N2λ+N1) dλ dτ = 0, (10) which can be arranged in the form of( dλ dτ )−1 = Q1(λ) Q2(λ) − τ λ , (11) where Q1 = ( 3λ2 + 2M2λ+M1 ) eλτ + 2N2λ+N1 Q2 = λ ( λ3 +M2λ 2 +M1λ+M0 ) eλτ Taking real component of ( dλ dτ )−1 at τ = τ0, with λ = iw, we have Re [ dλ dτ ]−1 τ=τ0 = BRNR +BINI N2 R +N2 I , where BR = 3w2cosτ0w0 − 2M2w 2cosτ0w0 −M1sinτ0w0 − 2N2w 2 − τw4cosτ0w0 −M2w 2cosτ0w0 +M2w 3sinτ0w0 −M0wsinτ0w0; BI = 3w3cosτ0w0 + 2M2w 2sinτ0w0 −M1wcosτ0w0 −N1w + τw4sinτ0w0 +M2w 3cosτ0w0 +M2 wsinτ0w0 −M0wcosτ0w0; NR = −w5sinτ0w0 +M2w 4cosτ0w0 +M1w 3sinτ0w0 −M0cosτ0w0 NI = w5cosτ0w0 +M2w 4sinτ0w0 −M1w 3cosτ0w0 −M0w 2sinτ0w0 Observe that if C3 : BRNR + BINI 6= 0 holds, then Re [ dλ dτ ]−1 τ=τ0 6= 0. Following the workings above and the Hopf bifurcation theory in [2, 3, 5, 7], we have the theorem below Theorem 2. . If conditions C1 − C3 hold, then the endemic equilibrium E∗(S∗, V ∗, I∗) of the system (1) is locally asymptotically stable when τ ∈ [0, τ0]; the system undergoes a Hopf bifurcation at E∗(S∗, V ∗, I∗) when τ = τ0 and a family of periodic solutions bifurcate from E∗(S∗, V ∗, I∗). L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 847 4. Numerical simulations and Discussions In this section we have carried out the simulations to validate the analytical findings and illustrate the long term dynamics of system (1). The parameter values are the same as the ones used in [8] with only τ being varied with time as indicated in the figures. Figure (2) shows that the disease free equilibrium is globally asymptotically stable when Rv = 0.7692 which is clearly less than unity. From the figure, it can be clearly seen that I0 = 0. Figure (3) shows that the endemic equilibrium E∗(35.7321, 45.5913, 6.3217) Figure 2: Simulation of system (1) shows the global stability of the disease-free equilibrium when Rv = 0.7692 is locally asymptotically stable when τ ∈ [0, τ0 = 31.1725]. This is in line with Theorem 2 above. L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 848 Figure 3: The effects of ε on all classes with (a) ε = 0.0127, τ = 5.76 and Rv = 4.5672, (b) ε = 0.91855, τ = 1.257, RV = 1.7261 and (c) ε = 0.4123, τ = 3.1267, RV = 3.1672. In Figure (3)(a), we can clearly see that when ε = 0.0127, τ = 5.76 and Rv = 4.5672, the rate of infectives are quite high as compared to Figures (3)(b) and (c) when ε = 0.91855, τ = 1.257, RV = 1.7261 and ε = 0.4123, τ = 3.1267, RV = 3.1672 respectively. This is a proof enough that rotavirus infections can be easily contained by introducing very strong vaccines and reducing the vaccine delay time. L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 849 Figure 4: Time plots of S, V and I with τ = 36.125 > τ0 = 31.1725 Figure 5: Bifurcation diagrams of system (1) with respect to τ (a)S, (b) V and (c) I L. O. Onyango et al. / Eur. J. Pure Appl. Math, 13 (4) (2020), 840-851 850 Figure 4 that a family of periodic solutions bifurcate at E∗(35.7321, 45.5913, 6.3217). This phenomenon is also illustrated by Figure 5 5. Conclusion and Recommendation In this work, we have formulated a mathematical model for rotavirus incorporating time delay in vaccination. The disease free equilibrium has been proved to be globally stable. The endemic equilibria is proved to be locally stable whenever τ = 0 and undergoes a Hopf bifurcation if τ > 0. From the analytical and simulation results, we recommend that a strong vaccine with a short delay time should be introduced in order to effectively control rotavirus infections. As a future work, we propose that the stability and directions of Hopf bifurcations derived in this work should be established. 6. Declarations 6.1. Competing interests The authors declare that they have no competing interests. 6.2. Authorial Contribution The authors have contributed as follows: (i) Florence Adongo: Formulated the model (ii) Lawrence Onyango: analyzed the model (iii) Job Bonyo : Performed numerical Analysis (iv) : Ogada A. Elisha: performed numerical simulations (v) George Lawi: proof read the manuscript 6.3. Source of data The data used to perform numerical simulation is from Onyango et.al [8] Acknowledgements We the authors would like to sincere thank Dr. Patriciah Gathia for the support and encouragement she gave to us during the time we were busy doing the research. We can’t fail to appreciate the many lunches you bought to us during this time. May the Almighty bless you Doctor. REFERENCES 851 References [1] Rotavirus and symptoms. Center for Disease Control, sep 2016. [2] Omar Bazighifan and Clemente Cesarano. Some new oscillation criteria for second order neutral differential equations with delayed arguments. Mathematics, 7(7):619, 2019. [3] Omar Bazighifan and Clemente Cesarano. A philos-type oscillation criteria for fourth- order neutral differential equations. Symmetry, 12(3):379, 2020. [4] C Carlos-Chavez, F Zhilan, and W Huang. On the computation of and its role on global stability. Institute for Mathematics and its Application, 125, 2001. [5] Brian D Hassard, DB Hassard, Nicholas D Kazarinoff, Y-H Wan, and Y Wah Wan. Theory and applications of Hopf bifurcation, volume 41. CUP Archive, 1981. [6] Benjamin A Lopman, Virginia E Pitzer, Rajiv Sarkar, Beryl Gladstone, Manish Patel, John Glasser, Manoj Gambhir, Christina Atchison, Bryan T Grenfell, W John Ed- munds, et al. Understanding reduced rotavirus vaccine efficacy in low socio-economic settings. PloS one, 7(8):e41720, 2012. [7] Jerrold E Marsden and Marjorie McCracken. The Hopf bifurcation and its applica- tions, volume 19. Springer Science & Business Media, 2012. [8] Onyango Lawrence Omondi, Chuncheng Wang, Xiaoping Xue, and Owuor George Lawi. Modeling the effects of vaccination on rotavirus infection. Advances in Differ- ence Equations, 2015(1):381, 2015. [9] Virginia E Pitzer, Katherine E Atkins, Birgitte Freiesleben de Blasio, Thierry Van Ef- felterre, Christina J Atchison, John P Harris, Eunha Shim, Alison P Galvani, W John Edmunds, Cecile Viboud, et al. Direct and indirect effects of rotavirus vaccination: comparing predictions from transmission dynamic models. PloS one, 7(8):e42320, 2012. [10] Sine Leergaard Wiggers and Pauli Pedersen. Routh–hurwitz-liénard–chipart criteria. In Structural Stability and Vibration, pages 133–140. Springer, 2018.