Processing math: 100%
  Special Issues

Prediction of influenza peaks in Russian cities: Comparing the accuracy of two SEIR models

  • Received: 08 November 2016 Accepted: 06 April 2017 Published: 01 February 2018
  • MSC : 37N25, 65C20, 92C60

  • This paper is dedicated to the application of two types of SEIR models to the influenza outbreak peak prediction in Russian cities. The first one is a continuous SEIR model described by a system of ordinary differential equations. The second one is a discrete model formulated as a set of difference equations, which was used in the Baroyan-Rvachev modeling framework for the influenza outbreak prediction in the Soviet Union. The outbreak peak day and height predictions were performed by calibrating both models to varied-size samples of long-term data on ARI incidence in Moscow, Saint Petersburg, and Novosibirsk. The accuracy of the modeling predictions on incomplete data was compared with a number of other peak forecasting methods tested on the same dataset. The drawbacks of the described prediction approach and possible ways to overcome them are discussed.

    Citation: Vasiliy N. Leonenko, Sergey V. Ivanov. Prediction of influenza peaks in Russian cities: Comparing the accuracy of two SEIR models[J]. Mathematical Biosciences and Engineering, 2018, 15(1): 209-232. doi: 10.3934/mbe.2018009

    Related Papers:

    [1] K. E. Starkov, Svetlana Bunimovich-Mendrazitsky . Dynamical properties and tumor clearance conditions for a nine-dimensional model of bladder cancer immunotherapy. Mathematical Biosciences and Engineering, 2016, 13(5): 1059-1075. doi: 10.3934/mbe.2016030
    [2] Svetlana Bunimovich-Mendrazitsky, Yakov Goltser . Use of quasi-normal form to examine stability of tumor-free equilibrium in a mathematical model of bcg treatment of bladder cancer. Mathematical Biosciences and Engineering, 2011, 8(2): 529-547. doi: 10.3934/mbe.2011.8.529
    [3] Eugene Kashdan, Svetlana Bunimovich-Mendrazitsky . Hybrid discrete-continuous model of invasive bladder cancer. Mathematical Biosciences and Engineering, 2013, 10(3): 729-742. doi: 10.3934/mbe.2013.10.729
    [4] Guobing Lin, Cunming Zhang, Xuanyu Chen, Junwei Wang, Song Chen, Siyuan Tang, Tianqiang Yu . Identification of circulating miRNAs as novel prognostic biomarkers for bladder cancer. Mathematical Biosciences and Engineering, 2020, 17(1): 834-844. doi: 10.3934/mbe.2020044
    [5] Ji-Ming Wu, Wang-Ren Qiu, Zi Liu, Zhao-Chun Xu, Shou-Hua Zhang . Integrative approach for classifying male tumors based on DNA methylation 450K data. Mathematical Biosciences and Engineering, 2023, 20(11): 19133-19151. doi: 10.3934/mbe.2023845
    [6] OPhir Nave, Shlomo Hareli, Miriam Elbaz, Itzhak Hayim Iluz, Svetlana Bunimovich-Mendrazitsky . BCG and IL − 2 model for bladder cancer treatment with fast and slow dynamics based on SPVF method—stability analysis. Mathematical Biosciences and Engineering, 2019, 16(5): 5346-5379. doi: 10.3934/mbe.2019267
    [7] Tongmeng Jiang, Pan Jin, Guoxiu Huang, Shi-Cheng Li . The function of guanylate binding protein 3 (GBP3) in human cancers by pan-cancer bioinformatics. Mathematical Biosciences and Engineering, 2023, 20(5): 9511-9529. doi: 10.3934/mbe.2023418
    [8] Liang-Sian Lin, Susan C Hu, Yao-San Lin, Der-Chiang Li, Liang-Ren Siao . A new approach to generating virtual samples to enhance classification accuracy with small data—a case of bladder cancer. Mathematical Biosciences and Engineering, 2022, 19(6): 6204-6233. doi: 10.3934/mbe.2022290
    [9] Zekun Xin, Yang Li, Lingyin Meng, Lijun Dong, Jing Ren, Jianlong Men . Elevated expression of the MYB proto-oncogene like 2 (MYBL2)-encoding gene as a prognostic and predictive biomarker in human cancers. Mathematical Biosciences and Engineering, 2022, 19(2): 1825-1842. doi: 10.3934/mbe.2022085
    [10] Mahmoud El-Morshedy, Zubair Ahmad, Elsayed tag-Eldin, Zahra Almaspoor, Mohamed S. Eliwa, Zahoor Iqbal . A new statistical approach for modeling the bladder cancer and leukemia patients data sets: Case studies in the medical sector. Mathematical Biosciences and Engineering, 2022, 19(10): 10474-10492. doi: 10.3934/mbe.2022490
  • This paper is dedicated to the application of two types of SEIR models to the influenza outbreak peak prediction in Russian cities. The first one is a continuous SEIR model described by a system of ordinary differential equations. The second one is a discrete model formulated as a set of difference equations, which was used in the Baroyan-Rvachev modeling framework for the influenza outbreak prediction in the Soviet Union. The outbreak peak day and height predictions were performed by calibrating both models to varied-size samples of long-term data on ARI incidence in Moscow, Saint Petersburg, and Novosibirsk. The accuracy of the modeling predictions on incomplete data was compared with a number of other peak forecasting methods tested on the same dataset. The drawbacks of the described prediction approach and possible ways to overcome them are discussed.


    1. Introduction

    Dengue is a vector-borne disease spread by the female mos-quitoes Aedes aegypti and Aedes albopictus. The mosquitoes Aedes aegypti were originated from Africa but have now been spread in tropical, subtropical and temperate regions of the world. The fast growth in human population, uncontrolled urbanization and inadequate waste management systems have led to abundance of mosquito breeding sites [6]. These are the main causes which bring this global distribution of mosquitoes and consequently spread of the mosquito-borne diseases. Before 1970, dengue was reported from only nine countries, now it has been widespread in the areas of North America, South America, Africa and Southeast Asia [7]. It is important to control dengue disease as more than 50 million people are affected every year [27]. Interestingly, only the female mosquitoes can transmit the virus. The female mosquitoes bite the human as they require blood for reproduction process. On biting the dengue infected person, these female mosquitoes become infected and can transmit the dengue virus to another person [7]. One of the ways to control dengue infection is to control mosquito population. Mosquito population can be chemically controlled by the use of insecticides and / or biological controlled by wolbochia [9,18,26]. The human awareness campaigns can also reduce the spread of infection.

    One of the efficient ways to combat the infection is the use of sterile insect technique (SIT). In SIT, sterile male mosquitoes are released near the mosquito breeding sites. Female mosquitoes will not be able to fertilize when mating with these male mosquitoes. In this way, the mosquito population as well as the spread of infection is controlled. The idea of SIT was first conceived by American entomologist, Dr Edward F. Knipling and was successfully implemented to control the spread of screwworm fly in Florida [14]. After that, the SIT technique has been used for many flying insects by various countries [15,16]. To eradicate dengue infection, the government of China and Brazil have also released the sterile male mosquitoes to combat the dengue infection [10,11].

    In literature, some mathematical models on control of vector-borne diseases are available [1,2,3,4,9,13,18,19,20,21,23,25]. Few mathematical models have also been formulated to control dengue infection by SIT [3,4,25]. In particular, Esteva and Yang proposed a mathematical model to control the mosquito density by releasing sterile male mosquitoes [3]. Optimal control analysis to control the Aedes aegypti female mosquitoes by SIT strategy has been performed by Thomˊe, Yang and Esteva [25]. A pulsed spatial discrete time model using SIT method has also been proposed by Evans and Bishop [4]. They concluded that increased number of sterile will have no further benefit above a threshold. The effectiveness is highly reduced by density dependent mortality of sterile insects. Hendron and Bonsall have developed the n patch dengue model with two control strategies namely vaccination and vector control by SIT method. The dynamics of two serotypes of dengue viruses with effect of cross immunity and antibody-dependent enhancement have been considered without explicitly incorporating the male/female vector (mosquito) dynamics [8].

    In this paper, a host-vector model has been proposed to control the primary dengue transmission using SIT strategy. In SIT, the male mosquitoes (Aedes aegypti) are sterilized in the laboratory and are released into the environment to compete with wild male mosquitoes for mating with female mosquitoes. In section 2, a network model for n patches has been formulated. The model analysis of an isolated patch is carried out in section 3. In section 4, the model analysis for n patch network model has been carried out to explore the behavior of system. Conclusions are given in the last section.


    2. Formulation of host-vector model

    Consider A be the class of mosquitoes in aquatic stage (eggs, larvae or pupae). Let U1 and U2 be the normal male mosquitoes and sterile male mosquitoes population respectively. Let the female mosquitoes are divided in four compartments: F1, F2,F3 and F4. The unmated and uninfected mosquitoes are in F1 compartment. Let F2 be the class of fertilized female mosquitoes and are uninfected. The female mosquitoes (F3) mated with sterile male mosquitoes are unfertilized and will not reproduce. The fertilized female mosquitoes need multiple blood meals to complete the reproduction cycle. The fertilized female mosquitoes who bite infected individuals become infected and transfer to F4 class and play crucial role in the spread of dengue infection. The remaining female fertilized mosquitoes bite susceptible/exposed/recovered human population before laying eggs and remain uninfected. No separate class is considered for them. As such both female mosquitoes F2 and F4 contribute to aquatic stage. Let the susceptible (S), exposed (E), infected (I) and recovered (R) be the host population in the four compartments. Consider the SEIR dynamics in the host population interacting with susceptible female mosquitoes F2 and infected female mosquitoes F4. The model also incorporates the interstate transition between different classes of male and female mosquitoes. No vertical transmission has been assumed for vector dynamics. Let there be a network of n patches. The host-vector dynamics in each patch is shown in the schematic diagram in Figure 1. Let the subscript i represents the respective variable/parameter in the ith patch; i=1, ...n. The human migration is considered between patches while no interpatch migration between mosquitoes is considered. The superscripts S, E, I and R are used to incorporate the differential ability of migration in various compartments (e.g. mSij represents migration in susceptible compartment). The following n patch metapopulation model is formulated according to the schematic diagram, Figure 1.

    Figure 1. The transfer diagram representing the system (1)-(11) dynamics. The births, deaths and migration are not included.
    dSidt=ωiβ1iSiF4iμiSi+nj=1,jimSijSjnj=1,jimSjiSi (1)
    dEidt=β1iSiF4i(ki+μi)Ei+nj=1,jimEijEjnj=1,jimEjiEi (2)
    dIidt=kiEi(αi+μi)Ii+nj=1,jimIijIjnj=1,jimIjiIi (3)
    dRidt=αiIiμiRi+nj=1,jimRijRjnj=1,jimRjiRi (4)
    dAidt=ϕi(1AiCi)(F2i+F4i)(γi+di)Ai (5)
    dF1idt=piγiAiβ2iF1iU1iU1i+U2iβ2iF1iU2iU1i+U2id1iF1i (6)
    dF2idt=β2iF1iU1iU1i+U2iβ3iIiF2id1iF2i (7)
    dF3idt=β2iF1iU2iU1i+U2id1iF3i (8)
    dF4idt=β3iIiF2id1iF4i (9)
    dU1idt=(1pi)γiAid1iU1i (10)
    dU2idt=ω1id1iU2i (11)

    All model parameters are defined in Table 1 and are assumed to be non-negative. The initial conditions associated with the above system are:

    Si(0)0,Ei(0)0,Ii(0)0,Ri(0)0,0Ai(0)Ci,F1i(0)0,F2i(0)0,F3i(0)0,F4i(0)0,U1i(0)0,U2i(0)0
    Table 1. Parameters of the Model.
    Parameters Description of parameters
    α Human recovery rate
    β1 Transmission rate of infection from female mosquitoes to human
    β2 Mosquitoes mating rate
    β3 Transmission rate of infection from human to female mosquito
    γ Transition rate from aquatic stage to adult mosquito
    μ Natural death rate of human
    ω Birth rate of human
    ω1 Constant recruitment rate of sterile male mosquito
    ϕ Recruitment rate for aquatic mosquito
    C Carrying capacity for aquatic/adult mosquito
    d Natural death rate of mosquito at aquatic state
    d1 Natural death rate of mosquito
    k Rate at which exposed human become infectious
    mij Migration rate from patch j to patch i
    p Proportion of female mosquito
     | Show Table
    DownLoad: CSV

    In the next section, the dynamics of the model (1)-(11) is discussed for an isolated patch (i.e mSij=mEij=mIij=mRij=0).


    3. Model analysis for an isolated patch

    In absence of interpatch migration, the patches are decoupled and isolated. The dynamics of each isolated patch is given below where the subscript i is suppressed:

    dSdt=ωβ1SF4μS (12)
    dEdt=β1SF4kEμE (13)
    dIdt=kEαIμI (14)
    dRdt=αIμR (15)
    dAdt=ϕ(1AC)(F2+F4)(γ+d)A (16)
    dF1dt=pγAβ2F1U1U1+U2β2F1U2U1+U2d1F1 (17)
    dF2dt=β2F1U1U1+U2β3IF2d1F2 (18)
    dF3dt=β2F1U2U1+U2d1F3 (19)
    dF4dt=β3IF2d1F4 (20)
    dU1dt=(1p)γAd1U1 (21)
    dU2dt=ω1d1U2 (22)

    Let X(t)=(S(t),E(t),I(t),R(t),A(t),F1(t),F2(t),F3(t),F4(t),U1(t),U2(t))R11+. The system (12)-(22) with non-zero initial conditions can be written in the following form:

    dXdt=Q(X(t),t),X(0)=X00Q(X(t),t)=(Q1(X,t),Q2(X,t),Q3(X,t)....,Q11(X,t)) (23)

    3.1. Positivity and boundedness of the solution

    Proposition 1. The positive cone Int(R11+) is invariant for the system (23).

    Proof. Observe that the boundaries of non-negative cone R11+ are invariant for the system (23) and Q(X(t),t) is smooth enough. Applying the existence and uniqueness theorem [22] for differential equations, the system (23) will possess the positive solution.

    Proposition 2. Solutions of system (12)-(22) are bounded in the domain Ω= {(S,E,I,R,A,F1,F2,F3,F4,U1,U2)R11+|S+E+I+Rωμ, A(t)C,F1+F2+F3+F4+U1+U2ω1+γCd1}.

    Proof. Let N(t)=S(t)+E(t)+I(t)+R(t) be the total host population at time t. Adding corresponding equations for host dynamics gives

    dNdt=ωμNlim suptN(t)ωμ.

    For vector dynamics, it is clear that A(t)C t. Let us prove it by contradiction.

    Let t0 be the smallest value of t such that A(t0)=C. Assume A(t)C for t[0,t0) and A(t)>C, for t(t0,T) where, T<. Accordingly, ˙A(t)>0 for t(0,t0) and ˙A(t)<0 for t[t0,T) from (16). By mean value theorem ˙A(t)>0 for t[t0,T) which is contradiction.

    Let M(t)=F1+F2+F3+F4+U1+U2, now adding equations from (17)-(22) gives

    dMdt=ω1+γA(t)d1(F1(t)+F2(t)+F3(t)+F4(t)+U1(t)+U2(t))dMdtω1+γCd1M(t)lim suptM(t)ω1+γCd1.

    Therefore, the system (23) is bounded.


    3.2. Equilibrium states

    The equilibrium states for the system (12)-(22) are obtained by solving

    Q(Xt),t)=0.

    The trivial disease-free state ˘P1 always exists and it is given as

    ˘S=ωμ,˘U2=ω1d1,˘E=0,˘I=0,˘R=0,˘A=0,˘F1=0,˘F2=0,˘F3=0,˘F4=0,˘U1=0.

    Note that, this state is without native mosquitoes. Only the susceptible human and sterile male mosquitoes survive.

    Let us define the non-dimensional number T, the basic offspring number, as

    T=β2ϕpγ(β2+d1)d1(γ+d). (24)

    It represents the average number of secondary female mosquitoes produced by single female mosquito, Esteva and Yang [3]. To maintain the mosquito population in environment, T should be greater than 1. The existence of another non-trivial disease-free equilibrium state is established for T>1.

    Proposition 3. The system admits two non-trivial disease-free states P+2 and P2 when

    T>1and ω1<W(=(T1)2Cγ(1p)4T). (25)

    Proof. Let us denote P±2=(ˉS,ˉE,ˉI,ˉR,ˉA,¯F1,¯F2,¯F3,¯F4,¯U1,¯U2). For the disease-free state of the system (12)-(22), substituting I=0 gives

    ˉS=ωμ,ˉE=0,ˉR=0,¯F4=0and ¯U2=ω1d1. (26)

    The other state variables ¯F1,¯F2,¯F3 and ¯U1 at equilibrium level are obtained in terms of ˉA:

    ¯F1=pˉAγ(β2+d1),¯F2=β2pˉAγˉU1(ˉU1+ˉU2)(β2+d1)d1,¯F3=β2pˉAγˉU2d1(β2+d1)(ˉU1+ˉU2),¯U1=ˉAγ(1p)d1 (27)

    ˉA is the root of following quadratic polynomial:

    s(ˉA)=B1ˉA2+B2ˉA+B3=0B1=β2pϕγ(β2+d1)(γ+d)d1C,B2=β2ϕpγ(β2+d1)d1(γ+d)+1,B3=ω1γ(1p)

    By writing the quadratic s(ˉA) in terms of T gives

    s(ˉA)=TCˉA2(T1)ˉA+ω1γ(1p)=0. (28)

    Since p<1, no positive root of equation (28) is admissible for T1. Clearly, the quadratic has real positive roots when the condition (25) is satisfied. The roots are given as

    ˉA±=(T1)C2T(1±1ω1W); W=(T1)2Cγ(1p)4T. (29)

    Accordingly, there exists two non-trivial disease-free equilibrium states P2 and P+2 corresponding to ˉA and ˉA+ respectively under condition (25).

    These two roots ˉA and ˉA+ collapse at ω1=W. This gives a critical value T=T where the quadratic has the unique solution:

    (ˉA=)ˉA=(T1)C2T (30)
    T=(1+2ω1Cγ(1p))[1+(11(1+2ω1Cγ(1p))2)]>1 (31)

    Therefore, the unique non-trivial disease-free equilibrium point P2 exists at T=T.

    Let us denote the endemic states by P±3 (ˆS,ˆE,ˆI,ˆR,ˆA,^F1,^F2,^F3,^F4,^U1,^U2). Now, the conditions for the existence of the endemic states P±3 are explored. The equilibrium level of state variables are obtained as ˆS=ω(β1^F4+μ),ˆE=β1ˆSˆF4(k+μ),ˆI=β1kˆSˆF4(k+μ)(α+μ),ˆR=αˆIμ, ^F1=pγˆA(β2+d1), ^F2=((γ+d)ˆACβ1+μϕ(CˆA))d1(k+μ)(α+μ)ϕ(CˆA)(β1β3kω+β1d1(k+μ)(α+μ)), ^F3=β2pγˆAω1d1(β2+d1)(^U1+^U2), ^F4=β1β3kω^F2μd1(α+μ)(k+μ)β1d1(α+μ)(k+μ), ^U1=(1p)γˆAd1,^U2=ω1d1.

    Here, ˆA is the root of same quadratic polynomial s(ˆA)=0 as given in (28). Accordingly, the condition (25) is necessary for the existence of both the non-trivial disease-free states P±2 as well as the endemic states P±3. Another condition ^F4>0 is required for existence of P±3 states.

    Remark 1. It is to be noted that for the non-trivial disease-free states P±2 and endemic states P±3, the equilibrium densities of following variables are found to be the same:

    ˉA=ˆA,ˉF1=ˆF1,ˉU1=ˆU1,ˉU2=ˆU2

    Remark 2. Simplifying the expression for ˆF2 gives

    ^F2=β2γpˆA^U1(^U1+^U2)(β2+d1)(β3ˆI+d1).

    Comparing it with ¯F2 yields,

    ˉF2>ˆF2 (32)

    Let us define

    R+0=β1β3kωˉF2+μd1(α+μ)(k+μ) and R0=β1β3kωˉF2μd1(α+μ)(k+μ). (33)

    Further, for positive ^F4+ and ^F4, the following conditions should be satisfied:

    β1β3kωˆF2+μd1(α+μ)(k+μ)>1 and β1β3kωˆF2μd1(α+μ)(k+μ)>1 (34)

    Using (32) gives the conditions for the existence of the endemic states P±3.

    R+0>1 and R0>1 (35)

    Remark 3. It can be seen from (29) that ˉA+>ˉA. Similarly, it can be easily proved that ˉF2+>ˉF2. Accordingly, it can be concluded that:

    R+0>R0

    Let us define R0 as

    R0=max(R+0,R0)(>1). (36)

    Accordingly, the following proposition is established for the existence of endemic states P±3:

    Proposition 4. When condition (25) is satisfied, the two positive endemic equilibrium states P3 and P+3 will exist provided the equation (36) is satisfied.

    Consider the following choice of parameters [3]:

    d=0.05,d1=0.0714,p=0.5,γ=0.075

    The existence of P±2 and P±3 states depend on the T and W. The parameters β2 and ϕ are considered to be critical for their existence. The curve T=1 in Figure 2 divides the β2ϕ plane into two regions. On the right of this curve (T>1) both the equilibrium states P±2 and P±3 may exist depending on the other conditions. In particular, the states P±2 and P±3 will not exist on the left of the curve T=1 (T<1).

    Figure 2. The existence of the state P2 with respect to parameters β2 and ϕ.

    3.3. The basic reproduction number

    The basic reproduction number has been computed by next generation approach [12]. It is defined as the average number of secondary infections produced by single infected individual in susceptible population. Using the Remark 3, it is given as

    R0=β1β3kωˉF2+(α+μ)d1μ(μ+k). (37)

    3.4. Stability of trivial disease-free state (P1)

    For local stability, the Jacobian matrix J[P] of the system (12)-(22) about a state P is given as:

    [β1F4μ0000000β1S00β1F4kμ000000β1S000kαμ0000000000αμ00000000000j5,50j5,70j5,9000000pγβ2d10000000β3F200j7,6d1β3I00j7,10j7,1100000j8,60d10j8,10j8,1100β3F2000β3I0d1000000(1p)γ0000d100000000000d1]
    j5,5=γd+ϕC(F2+F4),j7,6=β2U1U1+U2,j7,10=β2F1U2(U1+U2)2,j7,11=β2F1U1(U1+U2)2,j8,6=β2U2(U1+U2),j8,10=j7,10,j8,11=j7,11j5,7=ϕ(1AC)=j5,9

    The eleven eigenvalues of the Jacobian matrix about the trivial disease-free state P1 are obtained as

    αμ,μ(multiplicity 2),β2d1,d1(multiplicity 5),μk,γd1.

    Since all the eigenvalues have negative real part, the state P1 will always be locally asymptotically stable.

    Proposition 5. The locally asymptotically stable trivial disease-free state P1 is also globally stable for

    T<1. (38)

    Proof. Consider the positive definite function L(A,F1,F2,F4) for arbitrarily chosen positive constants C1,C2, C3 and C4:

    L(A,F1,F2,F3,F4)=C1A+C2F1+C3F2+C4F4

    Computing its time derivative ˙L(A,F1,F2,F4) along the trajectories of system (12)-(22),

    ˙L(A,F1,F2,F4)=C1˙A+C2˙F1+C3˙F2+C4˙F4˙L(A,F1,F2,F4)=C1(ϕ(1AC)(F2+F4)(γ+d)A)+C2(pγAβ2F1d1F1)+C3(β2F1U1(U1+U2)β3IF2d1F2)+C4(β3IF2d1F4)˙L(A,F1,F2,F4)F2(C3d1C1ϕ)F4(C4d1C1ϕ)A(C1(γ+d)C2pγ)F1(C2(β2+d1)C3β2)

    For the derivative of the function L(A,F1,F2,F4) to be negative definite

    d1C3>ϕC1;d1C4>ϕC1;(γ+d)C1>pγC2;(β2+d1)C2>β2C3

    or

    C3C1>ϕd1,C4C1>ϕd1,C1C2>pγ(γ+d),C2C3>β2(β2+d1)

    Now, choosing C3=C4=1 and using above inequalities give

    β2ϕpγ(β2+d1)d1(γ+d)(=T)<1.

    Accordingly, the function L(A,F1,F2,F4) is a Lyapunov function for the condition (38). Since {P1} is the largest invariant set that contains the subset in which ˙L=0 for A=0,F1=0,F2=0,F4=0. By applying LaSalle's invariance principle [17], locally asymptotically stable disease-free state P1 is also globally asymptotically stable under (38).


    3.5. Local stability of non-trivial disease-free states (P±2)

    Proposition 6. For the local stability of non-trivial disease-free states P±2, the following conditions should be satisfied:

    R0(=β1β3¯F2+kω(α+μ)μd1(μ+k))<1 (39)

    and

    k(ˉA)(=D4)=TˉA2Cω1γ(1p)>0 (40)

    Proof. For the local stability of P±2 states, the 11 eigenvalues of Jacobian matrix of the system (12)-(22) about the states P±2 are given as: μ (multiplicity 2), -d1 (multiplicity 2) and the remaining eigenvalues are the roots of cubic polynomial (q3(λ)) and fourth degree polynomial (q4(λ)), given as

    q3(λ)=λ3+B1λ2+B2λ+B3=0
    B1=α+2μ+d1+kB2=μ(μ+2d1)+(μ+d1)k+α(μ+d1+k)B3=(α+μ)d1(μ+k)β1β3ˉF2kˉS

    The fourth degree polynomial is given as

    q4(λ)=λ4+D1λ3+D2λ2+D3λ+D4=0.
    D1=β2+d+3d1+ˉF2ϕC+γ>0D2=(ˉF2ϕ+C(d+d1+γ))(β2+3d1)C+β2d1>0D3=d1((γ+d)(β2+2d1)+d1(β2+d1)+(γ+d)(2β2+3d1)d1(ˉU1+ˉU2)(TˉA2Cω1d1γ(1p)(2β2+3d1)))D4=TˉA2Cω1γ(1p)

    By Routh-Hurwitz criteria, all the roots of q3(λ) will be negative provided

    R0(=β1β3¯F2+kω(α+μ)μd1(μ+k))<1.

    Further, all the roots of fourth degree polynomial q4(λ) will have negative real part iff

    Di>0,i=1...4 (41)
    D1D2D3>D23+D21D4. (42)

    Note that D1 and D2 are positive and also D4 is positive provided the condition (40) is satisfied. This condition also ensures the positivity of D3. The condition (42) is simplified and found to be satisfied for D4>0.

    All eigenvalues have negative real part provided the conditions (39) and (40) are satisfied simultaneously and this completes the proof.

    Corollary 1. The non-trivial disease-free state P2 remains unstable always while the state P+2 is found to be locally asymptotically stable under the condition (39).

    Proof. For the stability of the two non-trivial disease-free states P±2 explicitly namely P2 and P+2, the sign of k(A) and k(A+) are critical. The k(ˉA) is defined in (40). For the point A where A+ and A collapse, we have AAA+. Using the expression for A as given in (30), k(A) is found to be zero. Accordingly, k(A) will be of negative sign giving instability of P2 and k(A+) will be of positive sign. Hence, P2 is always unstable while P+2 will be locally asymptotically stable for (39) to be satisfied.

    Remark 4. The condition (40) is always satisfied for P+2. Therefore its stability depends on the condition (39).

    Remark 5. When condition (39) is not satisfied, the non-trivial disease-free state P+2 will also become unstable. Consequently, the endemic states P+3 and P3 will start to exist by (36).

    The Figure 3 is two parameter bifurcation diagram with respect to β1 and β3 showing the region of stability of the states P±2 as well as the existence of the P±3 states. Consider the following choice of parameters with varying transmission rates β1 and β3:

    Figure 3. Bifurcation diagram for the stability of the state P2.
    ω=0.039,β2=0.5,k=0.1667,ϕ=0.5,μ=0.00004,d=0.1096,d1=0.0714,α=0.7,p=0.5,ω1=0.97,γ=0.7,C=500

    In Figure 3, the curve R0=1 bifurcates the parameter plane β1β3 into two regions. The non-trivial disease-free state P+2 is locally stable in the region R0<1. The existence of endemic states P+3 and P3 is possible in the region R0>1.


    3.6. Numerical simulations

    The local stability of P±3 states could not be achieved analytically due to large and complex expressions. To discuss the local stability of P±3 states, the numerical simulations have been performed for the data given in Table 2.

    Table 2. [24,3,5].
    Parameters Parameters values
    α 0.3
    β1 0.02
    β2 0.7
    β3 0.03
    γ 0.075
    ω 0.002
    ϕ 5
    μ 0.0000456
    C 450
    d 0.05
    d1 0.0714
    k 0.1667
    p 0.5
     | Show Table
    DownLoad: CSV

    Further, the ratio between the carrying capacity to the net sterile male mosquitoes i.e. (Cω1/d1) is assumed to be 0.8 [3]. It is found that for the data set given in Table 2, the basic offspring number T(=19.06388689)>1 and the threshold R0 =8.179178686(>1). Accordingly, the states P3+(0.76383733, 0.01181548, 0.00656461, 50.31770846,394.752323, 19.19005978, 46.75477253,141.25412, 0.128961,207.32790,624.64985) and P3(15.250814, 0.008414931, 0.00467529, 35.836020, 79.020, 3.8413960, 2.3417136, 35.314431, 0.00460008, 41.502141, 624.64985) exist by condition (36). For the stability of P3 state, numerical simulations have been performed for the four initial conditions Y1, Y2, Y3 and Y4 (defined in caption) in the neighborhood of the state P3.

    It is observed that the trajectories starting from the initial conditions Y1 and Y2 converge to the point Y which is the projection of the state P+3 on the hyperplane. Moreover, it is checked that all other state variables also tend to corresponding values of P+3. Thus, the solution trajectories converge to P+3. Similarly, the trajectories starting from the initial conditions Y3 and Y4 converge to the state P1. It shows the unstable behavior of P3 state. The phase plot in 3D hyperplane SIF4 has been drawn in Figure 4a. Further, for the stability of P+3 state, numerical simulations have been carried out for the four initial conditions (defined in caption) in the neighborhood of the state P+3.

    Figure 4. (a) A 3D phase plot showing the unstable behavior for the state P3 for the four initial conditions Y1=(14, 0.007, 0.001, 34, 75, 3.5, 2.5, 32, 0.004, 39,600), Y2=(19, 0.009, 0.003, 32, 70, 3.0, 3.0, 36, 0.009, 35,615), Y3=(20, 0.009, 0.08, 30, 78, 3.9, 2.0, 38, 0.09, 42,620) and Y4=(15, 0.006, 0.05, 37, 72, 3.2, 1.5, 30, 0.02, 45,630). (b) A 3D phase plot showing the stable behavior for the state P+3 for the four initial conditions Z1=(0.4, 0.013, 0.005, 40,400, 18, 50,140, 0.15,200,600), Z2=(1.8, 0.025, 0.009, 51,410, 21, 43,145, 0.09,210,630), Z3=(1.5, 0.022, 0.008, 48,398, 14, 45,132, 0.11,194,615) and Z4=(0.8, 0.018, 0.007, 45,390, 17, 48,138, 0.18,198,610).

    It is found that starting with the initial conditions Z1, Z2, Z3 and Z4, the solution converges to the P+3 state. The phase plot in 3D hyperplane SIF4 has been drawn in Figure 4b.

    A bifurcation diagram has been drawn for the system (12)-(22) with respect to sterile male mosquitoes rate ω1 in Figure 5. The transmission rates from vector to host(β1) and host to vector (β3) are considered to be 0.002 nd 0.003 and rest of the parameter values are taken from Table 2. It is observed that the non-trivial disease-free states P+2 and P2 exist from ω1=0 to W=80.0699 and collapse in P2 at ω1=W(=80.0699). The expression for W is given in (29). For W>80.0699, only the state P1 exists and stable. Further, it is found that at ω1=ω1 =30, R0= 1.002868566. The value of ω1 can be obtained using the condition (37). Therefore, for ω1>30, the non-trivial disease-free state P+2 is found to be stable. However, for ω1<30, the endemic states P+3 and P3 start to exist while the disease-free state P+2 becomes unstable. Thus, for ω1<W, the control on natural mosquitoes will depend on the initial conditions and beyond this W, the control will work successfully irrespective of any initial condition.

    Figure 5. Bifurcation plot of system (12)-(22) with respect to ω1.

    3.7. Discussion

    In this paper, a host-vector model has been formulated to analyze the effect of SIT control. It is observed from the analysis that when basic offspring number (T) is less than one then all the natural mosquitoes will be eliminated and only the sterile mosquitoes will survive in the long run. On the other hand if T is greater than T (given in (31)), the two non-trivial disease-free states start to exist in which bigger one is found to be stable for R0<1 while the lower one is unstable always. At T=T, the non-trivial disease-free states collapse. However, for R0>1, two endemic states exist. Numerical simulations have been carried out to analyze the stability of endemic states. It is found that one state gets stable while other remains unstable for choice of relevant data from literature. To control the disease, a threshold has been computed for the rate ω1 at which the sterile mosquitoes be introduced. If the initial population size is in neighborhood of disease-free state then SIT control will work effectively and disease will die out for ω1<ω1<W. However, when ω1<ω1, the disease will persist, as after this value endemic states start to exist (shown in Figure 5).

    To study the effect of human migration in disease transmission, the model (1)-(11) has been analyzed in next section.


    4. n patch network model analysis


    4.1. Basic reproduction number for n patch model

    When no movement of human is considered, R0 is given as

    R0=max(R0i);i=1,2,3...nwhere,R0i=kiβ1iβ3iωiF20i(μi+αi)(ki+μi)μid1i

    When the patches are coupled due to movement of human population, the basic reproduction number for the network model (1)-(11) is computed below:

    Let us denote the susceptible human population and fertilized female population at disease free state for the ith node be S0i and F20i and define

    qi=β1iS0i and li=β3iF20i

    Further, define n×n diagonal matrices Dq, Dl, Dk and D1d as

    Dq= diag(q1,q2,....qn), Dl= diag(l1,l2,....ln), Dk= diag(k1,k2,....kn) and D1d= diag(d11,d12,....d1n)

    The matrices QE and QI are also defined as

    QE=[μ1+k1+1=nj=1mEj1mE12...mE1nmE21μ2+k2+1=nj=1mEj2...mE2n.................mEn1.......μn+kn+1=nj=1mEjn]

    and

    QI=[μ1+α1+1=nj=1mIj1mI12...mI1nmI21μ2+α2+1=nj=1mIj2...mI2n...................mIn1......μn+αn+1=nj=1mIjn]

    The basic reproduction number is computed by next generation approach [12]. The Jacobian matrices of the system (1)-(11) for the new infections (F) and transfer from one compartment to another (Y) are given below:

    F=[00Dq0000Dl0]

    and

    Y=[QE00DkQI000Dd1]
    FY1=[00DqD1d1000Q1EQ1IDkDlDlQ1E0]

    The dominant eigenvalue of the next generation matrix FY1 is the basic reproduction number for n patch model, Rn0(say)

    Rn0=ρ(Q1EQ1IDkDlDqD1d1). (43)

    Particularly, for single node (n=1), R0 is obtained as

    R0=k1β11β31S01F201(μ+α1)(k1+μ)d1.

    4.2. Numerical examples

    Let us perform the numerical experiments to control the disease in a network using SIT. For n=2 and n=3, the network topologies are given in Figure 6 and Figure 7 respectively. Let us assume that mSij=mEij=mIij=mRij=mij and ω1i and mij are varying in different patches. For simplicity, rest of the data is considered to be the same in all patches given in Table 3.

    Figure 6. The network topology for n=2.
    Figure 7. The network topology for n=3.
    Table 3. [24,3,5].
    Parameters Parameters values
    αi 0.5
    β1i 0.001
    β2i 0.7
    β3i 0.001
    γi 0.075
    ωi 0.029
    ϕi 5
    μi 0.0000456
    Ci 450
    di 0.05
    d1i 0.0714
    ki 0.1667
    pi 0.5
     | Show Table
    DownLoad: CSV
    Network with n=2

    Observe that R0=3.6819(>1) in all the isolated patches for the data set given in Table 3. Accordingly, the disease is endemic in patches. Let us apply SIT in patch-1 so that the disease is controlled in it. The role of network coupling in controlling the disease in the network is explored for different combinations of migration parameters mij in Table 4. Further, if SIT control is reduced in patch-1 i.e. ω11=60 and no control is applied in patch-2 then for different combinations of migration parameters, the cases have been discussed in Table 5. It may be noted from Table 4, case (a) that the disease is controlled in the first patch where SIT is applied. The migration from second patch may bring infection in first patch. If the adequate number of sterile mosquitoes are present in the first patch, this additional infection will be eliminated. Consequently, the network is disease-free and SIT is successful. When ω11 is reduced (Table 5), the number of sterile mosquitoes present in the first patch are not sufficient enough to control the disease. Accordingly, the disease persists in the network (case (e)). However, changing the coupling (migration) suitably may again control the disease in the network (case (f)).

    Table 4. When SIT is applied only in the patch-1 i.e.(ω11=70, ω12=0).
    Cases Migration in patch-1 Migration in patch-2 Conclusions
    (a) m12=0 m21=0 Isolated patches with R0 of patch-1 is 0.2515<1 (37). Disease is controlled in patch-1 only.
    (b) m12=0.005 m21=0.0006 R20 for two patch network model is 0.7678<1 (43). Two patch network model may be disease-free. The time-series confirms that the infection level tends to zero in both the patches Figure 8a.
    (c) m12=0.009 m21=0.0025 R20 for two-patch network model is 1.4763>1 (43). SIT method will not be able to control disease with this migration combination Figure 8b.
     | Show Table
    DownLoad: CSV
    Figure 8. (a) Time series for I1 (black colour) and I2 (grey colour) converge to zero. (b) Time series for I1 and I2 converge to endemic state. (c) Time series for I1 and I2 converge to endemic state. (d) Time series for I1 and I2 converge to disease-free state. (e)Time series for I1 (black colour), I2 (dotted line) and I3 (grey colour) for different patches converge to disease-free state. (f)Time series for I1, I2 and I3 converge to endemic state.
    Table 5. When SIT is applied only in the patch-1 i.e.(ω11=60, ω12=0).
    Cases Migration in patch-1 Migration in patch-2 Conclusions
    (d) m12=0 m21=0 Isolated patches with R0 of patch-1 is 0.2515<1. Disease is controlled in patch-1 only.
    (e) m12=0.005 m21=0.001 R20 for two-patch network model is 1.1778>1. SIT method will not be able to control disease with this migration combination. Disease may persist in the network Figure 8c.
    (f) m12=0.009 m21=0.001 R20 for two-patch network model is 0.7273<1. The network may now be disease free Figure 8d.
     | Show Table
    DownLoad: CSV
    Networkwithn=3

    In this case, the following two cases have been considered in Table 6: From the Table 6, it is observed that the basic reproduction number for case (g) is less than one. The disease may be controlled in 3 node network model for this case. By numerically solving the network model (1)-(11) for n=3, it is found that the infection level converges to zero as shown in Figure 8e. Further, for case (h), (Rn0)>1. The time series for infective population (I1, I2 and I3) have been drawn in Figure 8f. It is evident that the infection persists. Thus, disease is controlled in case(g), but there is a failure of SIT in case (h).

    Table 6. When SIT is applied only in the patch-1 i.e.(ω11=70, ω12=0 and ω13=0).
    Cases Migration in patch-1 Migration in patch-2 Migration in patch-3 Rn0 for network model
    (g) m12=0.0001, m13=0.002 m21=0.0001 m31=0.0001 0.6762(<1)
    (h) m12=0.0002, m13=0.0002 m21=0.001 m31=0.001 2.6813(>1)
     | Show Table
    DownLoad: CSV

    Keeping the results of above numerical experiments, the following conclusions can be drawn:

    It is not necessary to apply SIT in the whole network. The disease can be controlled in network by applying SIT in one patch only. This is possible with suitable coupling of patches. Further, this selective way of applying SIT is more economical (cost effective). The success of SIT depends on the coupling strength of the network (migration parameters mij) and the recruitment rate of sterile mosquitoes ω11. Accordingly, they are suggested to be critical parameters.


    5. Conclusion

    In this paper, n patch network model has been formulated to control disease transmission by using SIT. First, the dynamics of an isolated patch has been analyzed. The basic reproduction number has been computed. For the existence and stability of disease-free and endemic states, the two critical parameters namely basic offspring number (T) and basic reproduction number (R0) have been identified. The bifurcation diagram has been plotted to show the existence and stability regions of disease-free and endemic states for an isolated patch. The critical level of sterile male mosquitoes has been obtained to control the disease. For n patch model, the basic reproduction number has been computed. The numerical simulations have been performed by considering two and three nodes with different combinations of migration and SIT recruitment to study the effects in disease elimination and persistence. It is concluded that the cost effective way to control the disease in suitably coupled network is to apply SIT in one patch only.


    Acknowledgments

    We would like to thank IFCAM for visiting support for this research collaboration. Authors also would like to thank anonymous reviewers for their valuable comments. Their comments have enhanced the clarity and the quality of results and text.


    [1] [ O. Baroyan, U. Basilevsky, V. Ermakov, K. Frank, L. Rvachev and V. Shashkov, Computer modelling of influenza epidemics for large-scale systems of cities and territories, in Proc. WHO Symposium on Quantitative Epidemiology, Moscow, 1970.
    [2] [ R. Burger,G. Chowell,P. Mulet,L. Villada, Modelling the spatial-temporal progression of the 2009 A/H1N1 influenza pandemic in Chile, Mathematical Biosciences and Engineering, 13 (2016): 43-65.
    [3] [ CDC, Influenza signs and symptoms and the role of laboratory diagnostics, [online], http://www.cdc.gov/flu/professionals/diagnosis/labrolesprocedures.htm.
    [4] [ CDC, People with heart disease and those who have had a stroke are at high risk of developing complications from influenza (the flu), [online], http://www.cdc.gov/flu/heartdisease/.
    [5] [ J. -P. Chretien, D. George, J. Shaman, R. A. Chitale and F. E. McKenzie, Influenza forecasting in human populations: A scoping review, PloS one, 9 (2014), e94130.
    [6] [ A. D. Cliff, P. Haggett and J. K. Ord, Spatial Aspects of Influenza Epidemics, Routledge, 1986.
    [7] [ V. Colizza, A. Barrat, M. Barthelemy, A. -J. Valleron and A. Vespignani, Modeling the worldwide spread of pandemic influenza: Baseline case and containment interventions, PLoS Med, 4 (2007), e13.
    [8] [ S. Cook, C. Conrad, A. L. Fowlkes and M. H. Mohebbi, Assessing google flu trends performance in the united states during the 2009 influenza virus a (h1n1) pandemic PloS one, 6 (2011), e23610.
    [9] [ N. Goeyvaerts, L. Willem, K. Van~Kerckhove, Y. Vandendijck, G. Hanquet, P. Beutels and N. Hens, Estimating dynamic transmission model parameters for seasonal influenza by fitting to age and season-specific influenza-like illness incidence Epidemics, 13 (2015), p1.
    [10] [ I. Hall,R. Gani,H. Hughes,S. Leach, Real-time epidemic forecasting for pandemic influenza, Epidemiology and Infection, 135 (2007): 372-385.
    [11] [ D. He, J. Dushoff, R. Eftimie and D. J. Earn, Patterns of spread of influenza A in Canada, Proceedings of the Royal Society of London B: Biological Sciences, 280 (2013), 20131174.
    [12] [ K. S. Hickmann,G. Fairchild,R. Priedhorsky,N. Generous,J. M. Hyman,A. Deshpande,S. Y. Del Valle, Forecasting the 2013-2014 influenza season using wikipedia, PLoS Comput Biol, 11 (2015): e1004239.
    [13] [ A. Hyder, D. L. Buckeridge and B. Leung, Predictive validation of an influenza spread model PloS one, 8 (2013), e65459.
    [14] [ F. Institute, Research Institute of Influenza website, [online], http://influenza.spb.ru/en/.
    [15] [ Y. G. Ivannikov and A. T. Ismagulov, Epidemiologiya Grippa (The Epidemiology of Influenza), Almaty, Kazakhstan, 1983, In Russian.
    [16] [ Y. Ivannikov,P. Ogarkov, An experience of mathematical computing forecasting of the influenza epidemics for big territory, Journal of Infectology, 4 (2012): 101-106.
    [17] [ V. N. Leonenko,S. V. Ivanov, Fitting the SEIR model of seasonal influenza outbreak to the incidence data for Russian cities, Russian Journal of Numerical Analysis and Mathematical Modelling, 31 (2016): 267-279.
    [18] [ V. N. Leonenko,S. V. Ivanov,Y. K. Novoselova, A computational approach to investigate patterns of acute respiratory illness dynamics in the regions with distinct seasonal climate transitions, Procedia Computer Science, 80 (2016): 2402-2412.
    [19] [ V. N. Leonenko,Y. K. Novoselova,K. M. Ong, Influenza outbreaks forecasting in Russian cities: Is Baroyan-Rvachev approach still applicable?, Procedia Computer Science, 101 (2016): 282-291.
    [20] [ V. N. Leonenko,N. V. Pertsev,M. Artzrouni, Using high performance algorithms for the hybrid simulation of disease dynamics on CPU and GPU, Procedia Computer Science, 51 (2015): 150-159.
    [21] [ D. C. Liu,J. Nocedal, On the limited memory bfgs method for large scale optimization, Mathematical programming, 45 (1989): 503-528.
    [22] [ A. Romanyukha,T. Sannikova,I. Drynov, The origin of acute respiratory epidemics, Herald of the Russian Academy of Sciences, 81 (2011): 31-34.
    [23] [ L. A. Rvachev,I. M. Longini, A mathematical model for the global spread of influenza, Mathematical Biosciences, 75 (1985): 1-22.
    [24] [ J. Shaman, V. E. Pitzer, C. Viboud, B. T. Grenfell and M. Lipsitch, Absolute humidity and the seasonal onset of influenza in the continental United States, PLoS Biol, 8 (2010), e1000316.
    [25] [ J. Tamerius,M. I. Nelson,S. Z. Zhou,C. Viboud,M. A. Miller,W. J. Alonso, Global influenza seasonality: Reconciling patterns across temperate and tropical regions, Environmental Health Perspectives, 119 (2011): 439-445.
    [26] [ J. Truscott,C. Fraser,S. Cauchemez,A. Meeyai,W. Hinsley,C. A. Donnelly,A. Ghani,N. Ferguson, Essential epidemiological mechanisms underpinning the transmission dynamics of seasonal influenza, Journal of The Royal Society Interface, 9 (2011): 304-312.
    [27] [ S. P. van Noort,R. Águas,S. Ballesteros,M. G. M. Gomes, The role of weather on the relation between influenza and influenza-like illness, Journal of Theoretical Biology, 298 (2012): 131-137.
    [28] [ C. Viboud,O. N. Bjornstad,D. L. Smith,L. Simonsen,M. A. Miller,B. T. Grenfell, Synchrony, waves, and spatial hierarchies in the spread of influenza, Science, 312 (2006): 447-451.
    [29] [ WHO, Influenza (seasonal). Fact sheet No. 211, March 2014. , [online], http://www.who.int/mediacentre/factsheets/fs211/en/.
    [30] [ WHO, Surveillance case definitions for ILI and SARI, [online], http://www.who.int/influenza/surveillance_monitoring/ili_sari_surveillance_case_definition/en/.
    [31] [ R. Yaari, G. Katriel, A. Huppert, J. Axelsen and L. Stone, Modelling seasonal influenza: The role of weather and punctuated antigenic drift, Journal of The Royal Society Interface, 10 (2013), 20130298.
    [32] [ W. Yang, B. J. Cowling, E. H. Lau and J. Shaman, Forecasting influenza epidemics in hong kong, PLoS Comput Biol, 11 (2015), e1004383.
  • This article has been cited by:

    1. Marom Yosef, Svetlana Bunimovich-Mendrazitsky, Mathematical model of MMC chemotherapy for non-invasive bladder cancer treatment, 2024, 14, 2234-943X, 10.3389/fonc.2024.1352065
    2. Christoph Sadée, Stefano Testa, Thomas Barba, Katherine Hartmann, Maximilian Schuessler, Alexander Thieme, George M Church, Ifeoma Okoye, Tina Hernandez-Boussard, Leroy Hood, Ilya Shmulevich, Ellen Kuhl, Olivier Gevaert, Medical digital twins: enabling precision medicine and medical artificial intelligence, 2025, 25897500, 100864, 10.1016/j.landig.2025.02.004
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4979) PDF downloads(717) Cited by(10)

Article outline

Figures and Tables

Figures(10)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog