Authors: Kazi Mehedi Mohammad, Taufiquar Khan, Md Kamrujjaman
Categories: Article, Lévy jump, Influenza, SDE, Seasonality, Random fluctuations
Source: Infectious Disease Modelling
Authors: Kazi Mehedi Mohammad, Taufiquar Khan, Md Kamrujjaman
Forecasting influenza outbreaks remains a significant challenge due to the complexity of disease transmission and the influence of environmental and behavioral factors. Traditional models based solely on the basic reproduction number (R0) often fall short in capturing the full scope of outbreak dynamics.
In this study, we employ a seasonally adjusted SEIRT model incorporating stochastic differential equations (SDEs), including Brownian motion and Lévy jump processes, to simulate random and abrupt fluctuations in transmission. A branching process approximation is used to evaluate the probability of an epidemic under the influence of seasonal variability and stochastic perturbations. The model is calibrated using weekly influenza case data from Mexico, with noise components estimated from publicly available CDC [1] and WHO [2] surveillance data.
Simulation results show that the inclusion of stochastic effects and periodic transmission rates significantly enhances the model's accuracy in reflecting real-world epidemic dynamics. Numerical comparisons between deterministic, Brownian-based, and Lévy-based scenarios reveal that both the initial state of the exposed or infectious subpopulation and the seasonal transmission patterns are critical to determining outbreak probabilities. Results indicate that seasonal transmission rates and stochastic effects significantly alter epidemic probabilities, with Lévy processes capturing abrupt outbreak dynamics more accurately than deterministic models.
The findings underscore that deterministic models may underestimate epidemic risk when they overlook random and sudden changes in contact rates or disease introduction. The proposed stochastic modeling framework yields a deeper understanding of influenza transmission dynamics by incorporating uncertainty and seasonal variability, thereby supporting more informed and effective public health decision-making.
A stochastic model's parameters or state variables are arbitrary, and the outcomes are unknown beforehand (Ang, 2007; Yang, 2016). Trans-disciplinary research is the most crucial component in enabling everyone to lead normal, healthy lives in contemporary society. Mathematical predictive modelling is one of the most crucial methods for researching and comprehending how epidemics behave. It helps decision-makers make important decisions and arrange their days ahead of time.The 2009 pandemic caused by the H1N1 virus received extensive media attention (Kamrujjaman & Mohammad, 2025; Kanyiri et al., 2018; Mohammad et al., 2025; Rosyada & Hariyanto, 2019). A CDC report states that the H1N1 virus is primarily spread by respiratory droplets emitted when an infected person coughs, sneezes, or speaks. Those that are in close contact to the virus may be exposed by inhaling these droplets, which could lead to the virus entering their respiratory system. Furthermore, the virus can be disseminated by touching the mouth, nose, or eye after coming into contact with an infected surface. The WHO claims that the virus may inadvertently spread since it can spread even when there are no symptoms. People infected with H1N1 are frequently contagious for up to seven days following infection onset and up to one day prior to the onset of symptoms. However, the duration of contagiousness varies from person to person, especially among unyielding individuals (Khanh, 2016; Krishnapriya et al., 2017; Leonenko & Ivanov, 2016).
Influenza has evolved into a complex parasitic disease on a global scale. In addition to spreading from person to person, the H1N1 virus can also spread from pigs to humans. Direct contact with sick pigs or contact with areas polluted by swine excrement or bodily fluids can also result in transmission. However, the main cause of H1N1 outbreaks is persistent human-to-human transmission. If you wish to stop the spread of H1N1, you must practise good respiratory hygiene. This means that when you cough or sneeze, cover your mouth and nose with an elbow or tissue and dispose of old tissues properly (Kamrujjaman et al., 2021; Lee & Chowell, 2017). It is crucial to often wash your hands with soap and water or use alcohol-containing hand sanitisers to reduce the risk of contamination. Staying at home if you have flu-like symptoms and avoiding close contact with sick people are important protective measures (Akter Akhi et al., 2023; Pitchaimani & Krishnapriya, 2016).
According to recent studies, a given group of diseases has a 1−1R0i chance of experiencing a severe disease epidemic, where i represents the infected persons and R0>1 is the threshold quantity. She was particularly well-known in the domains of continuous-time Markov chains and stochastic differential equations. It is increasingly vital to find leading indicators for respiratory diseases, which represent a global health burden, in order to prevent disease epidemics. Increasing human immunity can help lower the number of exposed and ill individuals (Khanh, 2014; Maji & Ghosh, 2024).
Seasonal changes or the passage of time might affect how infectious illnesses behave. The likelihood of a disease outbreak might therefore happen on a regular basis and be impacted by seasonal changes as well as the original population impacted. Estimates from the SDE, Lévy leap, and branching process approximations may not be accurate if there are less than 100 populations. With different parameter values, the CTMC model can be used to estimate the probability of an outbreak and the initial number of COVID-19 infections (Akter Akhi & Kamrujjaman, 2023; Gray et al., 2011; Mehedi Mohammad et al., 2024; Niu et al., 2021).
A good strategy for preventing H1N1 infection is receiving the right care. The seasonal influenza vaccination should be administered annually to all people, but it should be given more frequently to individuals who are more susceptible to problems, such as young children, the elderly, pregnant women, and people with underlying medical disorders. The seasonal influenza vaccination frequently offers protection against H1N1 in addition to the virus. The fact that the information provided here is based on knowledge that is up to date as of September 2021 is more significant. Reputable health organisations and official government websites are the best places to find current and reliable information about influenza H1N1 (Alcaraz & Vargas-De-León, 2012; Modnak & Wang, 2017).
While existing literature (Ang, 2007) extensively discusses the use of a simple stochastic differential equation in the modeling of an epidemic along with real data. Our study aims to address this gap by investigating Branching process approximation with infinitesimal transition probabilities by Kolmogorov differential equations. Most studies (Allen, 2017; Allen & Lahodny Jr, 2012; Gray et al., 2011; Yang, 2016) focus on simple formulation, difference equation, Numerical simulation. However, there is limited analysis exploring the treatment implications, particularly in different waves. Our study endeavors to bridge this gap by examining the Influenza model with seasonality effect and probability of disease outbreak. We have also carried out qualitative study of stochastic model with treatment effect based on effective reproduction number R0. While numerous studies (Arruda et al., 2021; Bayram et al., 2018; Rao et al., 2012; Wang et al., 2017) have investigated the impact parameters by implementing average value, there is limited research exploring its effectiveness seasonal waves, such as sinusoidal curve. Our study seeks to fill this gap by examining seasonal variation with sensitive parameters, also the efficacy normally distributed parameters, contributing to a more comprehensive understanding to mitigate the pick level of infection cases. Existing literature (Elhiwi, 2021; Niu et al., 2021; Srivastav et al., 2022; Zhang et al., 2020) primarily focuses SDE solution with random diffusion. Meanwhile, we have compared the Lévy jump scenario with the standard SDE model to assess its impact on the SEIRT disease model. However, there is a lack of empirical research examining the noise term in each compartment based on density and parameter variation. Our study aims to address this gap by conducting a case study analysis in Mexico with regional basis real data to identify parameters effect in threshold quantity R0, meanwhile we also estimated noise term σ and σ1 bases in two different waves by data fitting with suggested model.
In order to account for the intrinsic stochasticity in illness transmission dynamics and to recognize the random fluctuations that impact the virus's propagation, influenza models contain Brownian motion and Lévy jump. By giving a more accurate representation of the uncertainties and unpredictabilities in the transmission process, this enhances the model's ability to replicate the complex and dynamic nature of influenza epidemics.
Seasonality effects recognize the impact of environmental factors like climate and human behavior on the propagation of the virus and are incorporated into influenza models to account for seasonal variations in illness transmission. By incorporating seasonality, influenza models become more realistic and more accurately depict observed patterns of illness transmission at different times of the year, which is beneficial for public health planning and intervention tactics.
We have studied modified five-compartment stochastic and deterministic models, including the mathematical model known as Susceptible-Exposed-Infectious-Treatment-Removal (SEIRT). This model is used to predict future basic reproduction number and infection cases as well as recent trends, recent probability outbreaks, seasonality effect, persistence nature, and extinction of disease. The objectives of this paper •Incorporating Brownian motion, Lévy jump scenario into the influenza illness model allows for a more accurate portrayal of uncertainty and unpredictability in the virus's propagation by accounting for the intrinsic randomness in disease transmission dynamics.•Enhancing the model's faithfulness to known real-world patterns by including seasonality effects into the model to represent fluctuations in disease transmission rates throughout different periods of the year, taking into account extrinsic factors like climate and human behavior.•Developing a mathematical model that incorporates seasonality, Lévy jump, Brownian motion to adapt to dynamic environmental conditions and provide a more thorough knowledge of how the influenza virus spreads over time in response to changing parameters.•Combining seasonality, Lévy jump and Brownian motion effects to improve the SEIRT model's predictive accuracy. This would guarantee that the model's predictions closely match empirical data and offer more trustworthy information for public health planning.•Offering a more thorough and comprehensive representation of the spread of disease by enabling the SEIRT model to incorporate both random fluctuations (Brownian motion) and temporal changes (seasonality) to capture the intricate dynamics of influenza transmission.
The incorporation of Brownian motion and seasonality into the SEIRT influenza model significantly enhances its ability to replicate real-world disease dynamics. These additions effectively capture both the seasonally variable nature of disease transmission and the inherent randomness associated with the spread of infection. By integrating Brownian motion, Lévy jumps, and seasonality, the model achieves a more accurate estimation of influenza transmission. This leads to an improved fit with actual epidemiological data, making the model a powerful and reliable tool for forecasting the future course of the disease. Moreover, the inclusion of these stochastic elements and seasonal variations renders the model more adaptable to changing environmental conditions. This adaptability is critical for understanding how the influenza virus responds to shifting external factors over time. The model also captures intricate features of influenza dynamics by accounting for both regular seasonal patterns and random fluctuations. This detailed representation provides a deeper and more comprehensive understanding of how the virus spreads across different environments and populations. A comparative analysis of the effects of Brownian motion and Lévy jumps within the SEIRT model reveals the significant role stochastic influences play in the transmission of influenza. Gaining insight into how random events impact transmission patterns is essential for informing and shaping future public health strategies. Finally, the enriched model, incorporating seasonality, Lévy jumps, and Brownian motion, offers critical insights that support the refinement of intervention strategies. With enhanced predictive precision, public health authorities can design more targeted and effective measures to control and mitigate the impact of influenza outbreaks and pandemics.
The structure of this paper is as in Section 2, we define the basic mathematical model and its formulation with compartmental flow diagram. The fixed points, the basic reproduction number and stability analysis are presented in Appendix A. Then in Section 3, formulation of Brownian motion and numerical results including noise in parameters are presented elaborately. In Section 4, several qualitative behaviour of the model, including existence and uniqueness of positive solution, stochastic disease free dynamics, stochastic endemic dynamics; extinction, persistence of disease by nature of R0, impact of treatment rate to reduce the outbreak is analyzed theoretically and graphically. In Section 5.1, the ODE model and the fundamental reproduction number are analyzed in relation to influenza spread in both seasonal and nonseasonal environments. For the differential equations, a time-nonhomogeneous stochastic process is also discussed in Section 5.2. The branching process approximation and numerical illustration are investigated in Sections 5.3, 5.4. This illustrates the seasonality of parameters in several compartments using normally distributed characteristics. Further, in Section 6 formulation of model with Lévy jump effects along with stability analysis is presented. Basen on this, the comparison scenario of Lévy jump and Brownian Motion is also investigated. Moreover, in Section 7 a case study in Mexico with comaprison to our SDE model is depicted. Finally, a brief and practical conclusion is offered in Section 8.
In 2017, Chairat presents a mathematical model of influenza transmission that serves as the foundation for many other infectious illness models. Chairat regarded affected and unaffected individuals as members of a fixed-size human population, assuming an exponentially distributed infectious time (Modnak & Wang, 2017). After that, Alcaraz modified this model, which we took into account for our further influenza research. It is assumed that the human populations in the Chairat and Alcaraz models are evenly distributed (Alcaraz & Vargas-De-León, 2012).
The entire human population is separated into five epidemiological susceptible (S), exposed (E), infected (I), recovered (R), and treated (T) individuals. This division is made without altering the system's generality or dynamical behaviour. That means modification of conventional SEIR epidemic process is applied to simulate for human infection (Martcheva, 2015). Let us define Y=(Y1,Y2,Y3,Y4,Y5)=(S,E,I,R,T). Then the mathematical model is given in the compact form as below,(2.1)Y′=F(Y,t),Y(0)=Y0.where, the initial conditions, Y0=(S0,E0,I0,R0,T0), and,F(Y,t)=[Λ−(β1E+β2I)S−μS,(β1E+β2I)S−(α+μ)E,αE−(μ+δ+γ+γ1)I,γI−μR,γ1I−μT].where, the total population for the model is defined byN(t)≡S(t)+E(t)+I(t)+R(t)+T(t).We have examined a new strain of influenza A belonging to the H1N1 subtype in our investigation. Based on two basic characteristics of the viral family, influenza virus, an enveloped virus belonging to the Orthomyxoviridae family, has a unique capacity for genetic diversity (Alcaraz & Vargas-De-León, 2012; Martcheva, 2015). In addition to coughing and sneezing, this virus spreads from sick to healthy people through contact. Subsidy populations are exposed following the latent period, and exposed compartments develop into infectious compartments. In order to control the illness outbreak, some members of the infected class then recover and go on to the treatment class.
Here, Λ is the population's recruitment/birth rate. For the hosts, μ and δ stand for the number of natural and disease-induced death rates per unit of time, respectively. γ is the disease recovery rate. Assume that the vulnerable person comes into contact with the exposed and contaminated areas. The transmission rate of infection from exposed to susceptible is β1, while the rate at which healthy humans contract the infection through contact with an infected class is β2. As a result, the infection rate for healthy hosts from mixing with the E and I compartments is β1ES and β2IS, respectively. The rate of infection progression from exposed to infected class is shown by α. Additionally, the rate of receiving treatment from the infected class is represented by γ1. Here, we took into account time per unit on a weekly basis. All of the model parameters are shown in Table 1. The flow diagram is shown in Fig. 1.Table 1Model parameters and state variables together with their descriptions.Table 1NotationDefinitionΛRate of recruitment/birth in S classβ1The likelihood of disease transmission with S and E**β2The likelihood of disease transmission with S and I**αRate of progression of infection from E to I class.γRecovery rate.γ1Rate of treatment of I class to the treated class T.μNatural death rate.δDisease induced death rate in I compartment.Fig. 1Diagram of the influenza SEIRT model's compartments.Fig. 1
The fact that the restored population is no longer sensitive to the decreased SEIRT model was taken into account in our investigation. Therefore, if the patients are able to recover from their sickness, it will be permanent. Consequently, we have decided to consider models of the SEIRT type.
When considering an epidemic model, Brownian motion refers to the stochastic or random fluctuations observed in the dynamics of the disease spread. When applied to an epidemic model, Brownian motion introduces randomness or uncertainty into the system, acknowledging that real-world epidemics are subject to numerous unpredictable factors. These factors can include individual variations in contact patterns, variations in disease transmission rates, or the occurrence of random events that affect the spread of the disease (Allen, 2017; Allen & Lahodny Jr, 2012). Incorporating Brownian motion into an epidemic model allows for the simulation of more realistic and dynamic scenarios. It helps capture the inherent variability and unpredictability of disease transmission, which can be particularly important when studying the long-term behaviour of an epidemic or making predictions about future trends (Ang, 2007; Gray et al., 2011; Mahmud et al., 2022).
By including random noise through Brownian motion, the model can account for the stochastic nature of disease dynamics, enabling researchers to assess the impact of uncertainties and better understand the range of possible outcomes. This can support the formulation of policies, the distribution of resources, and public health intervention decision-making processes during epidemics (Bayram et al., 2018). To incorporate Brownian motion with random noise mathematically in the epidemic model, we can introduce stochasticity to the differential equations by adding random variables. Brownian motion can be represented as a stochastic differential equation (SDE) and can be discritized using Euler-Maruyama method. Here is the procedure how we have applied Brownian motion by mean variance with random noise mathematically to the SEIRT model (1) as displayed in Section 2:(1)Define the deterministic model equations which is given to (1).(2)Introduce random noise to the equations by adding Browninan motion term to each equation. Let's denote the random variables as dW1, dW2, dW3, dW4 and dW5 which represent the Brownian motion increments associated with each equation.dS(t)=Λ−(β1E(t)+β2I(t))S(t)−μS(t)dt+σ1dW1(t),dE(t)=(β1E(t)+β2I(t))S(t)−(α+μ)E(t)dt+σ2dW2(t),dI(t)=αE(t)−(μ+δ+γ+γ1)I(t)dt+σ3dW3(t),dR(t)=γI(t)−μR(t)dt+σ4dW4(t),dT(t)=γ1I(t)−μT(t)dt+σ5dW5(t),(3)Specify the parameters which is given in Table 1 and their values given in Table 3. Let's denote the stochastic terms by σ1 = σS(t), σ1 = σE(t), σ3 = σI(t), σ4 = σR(t) and σ5 = σT(t) which are the variance of parameters representing the intensity of the Brownian motion for each equation of the SEIRT model.(4)Discritize the equations using the Euler-Maruyama method. Assuming a time step size of Δt, the discretized equation S(t+Δt)=S(t)+[Λ−(β1E(t)+β2I(t))S(t)−μS(t)]Δt+σ1ΔtdW1E(t+Δt)=E(t)+[(β1E(t)+β2I(t))S(t)−(α+μ)E(t)]Δt+σ2ΔtdW2I(t+Δt)=I(t)+[αE(t)−(μ+δ+γ+γ1)I(t)]Δt+σ3ΔtdW3R(t+Δt)=R(t)+[γI(t)−μR(t)]Δt+σ4ΔtdW4T(t+Δt)=T(t)+[γ1I(t)−μT(t)]Δt+σ5ΔtdW5where, ΔW1, ΔW2, ΔW3, ΔW4 and ΔW5 are independent random variables following a standard normal distribution (N(0, Δt)). More, importantly, the discretized equations incorporate the random noise through the increments.
ΔW1, ΔW2, ΔW3, ΔW4 and ΔW5.
In these equations, the terms σS(t), σE(t), σI(t), σR(t) and σT(t) represent the stochastic noise added to each respective compartment. These terms follow a Brownian Motion process with mean having zero and variance denoted by D, where D represents the diffusion coefficient. The assumptions made and the particular features of the system would determine the precise type and size of noise terms. The stochastic terms σS(t), σE(t), σI(t), σR(t) and σT(t) can be simulated as random variables at time step using appropriate methods like Monte-Carlo simulations on stochastic differential equations.
It is noted that, incorporating stochastic terms into the SEIRT model introduces randomness into the system, allowing for the modeling of uncertainty and variability (Allen, 2017; Allen & Lahodny Jr, 2012). For this thesis work, the choice of diffusion coefficient and the modeling of the noise terms required careful consideration and may involve calibration with emperical data or expert knowledge.
In this section, we have presented the scenario of each compartments of our model (1) graphically, by applying Browinan motion with mean-increment = 0.5, standard deviation = 5. The initial conditions are supposed as
S(0) = 990, E(0) = 5, I(0) = 4, R(0) = 1, T(0) = 0. In Fig. 2(a), we visualize that when more people are exposed to the disease, the number of susceptible people gradually declines. Along with time span from week 0 to 20, the density of susceptible compartmental population progressively decreases from 990 to 110. After that the noise occurs and susceptible population slightly progress to the level around 200 to 400. This reflects the scenario of reduction of infection rate. However, in the context of Brownian motion, the fluctuations in the susceptible compartment will be influenced by the mean and standard deviation parameters. The mean of 0.5 suggests a slight downward trend in the number of susceptible individuals over time, while the standard deviation of 5 indicates significant random fluctuations around the mean. In our stochastic framework, an increase in population density amplifies the intensity of random fluctuations in disease transmission. This imply that the healthy population evolves dynamically through recruitment, infection, and natural mortality.Fig. 2The curves depict five instances of brownian motion in influenza models, wherein (a) susceptible, infected, and recovered populations behave, and (b) exposed and treated populations behave. The parameter values are β1 = 0.35, β2 = 0.35, γ = 0.08, γ1 = 0.05, α = 0.25, μ = 0.02, δ = 0.01, mean-increment = 0.5, standard-deviation = 5.Fig. 2
Fig. 2(a) illustrates that the infected population increases initially and then declines as individuals recover and leave the infectious compartment. We notice that, in time interval 0 to 25 weeks, the infection density strikes to maximum level 400, corresponding to total human 1000. After the progress of time, the curves slightly started to decrease. In is noticed that, in weeks 58 to 70, infection rate lies around 0 level, that indicates extinction of disease from community. After week 80 the density of infected class fluctuates below 200. In the context of Brownian motion, the fluctuations in the infected compartment will be influenced by the mean and standard deviation parameters. The mean of 0.5 suggests a slight upward trend followed by a downward trend in the number of infected individuals over time, while the standard deviation of 5 indicates significant random fluctuations around the mean with respect to the density.
Fig. 2(a) depicts that, the recovered individuals increases over time as individuals move out of the infected compartment. In the context of Brownian motion, the fluctuations in the recovered compartment will be influenced by the mean and standard deviation parameters. The mean of 0.5 suggests a slight upward trend in the number of recovered individuals over time, while the standard deviation of 5 indicates significant random fluctuations around the mean. We observe that, along weeks 0 to 40 when infection gradually decrease, that results into swift increase in recovered compartment. The recovered population is at the 500–600 level after 40 weeks, which guarantees that the exposed and infected density is under control.
Moreover, Fig. 2(b) reveals that, with the progression of time span, gradual increase occurs in exposed class up to week 0 to 20. Exposed density reaches to its peak level 340 in week 21. After the vaccination and treatment rate increases, a swift decreasing scenario observed in exposed class. After 45 weeks, exposed population converges to the o level. However, the mean of 0.5 indicates a modest upward trend in the number of people exposed over time, while the standard deviation of 5 indicates significant random fluctuations around the mean.
By analysing Fig. 2(b) result, it indicates the rising tendency of treated compartment in time intervals 0 to 45. After weeks 60, the fluctuation level bounded between 320 and 400, while then total human is 1000. This indicates that, over a long time horizon, herd immunity effects lead to a large proportion of the population transitioning into the recovered and treated compartments, which contributes to the control of the disease burden. In our study, it is effective where a portion of infected individuals requires treatment. In the context of Brownian motion, the fluctuations in the treatment compartment will be influenced by the mean and standard deviation parameters which vary proportionally with the treated density.
In context of Brownian motion, the parameters can also vary with noise (Arruda et al., 2021; Gray et al., 2011; Niu et al., 2021). In this section we have presented some scenario in parameters with diffusion rate, σ = 0.1. Fig. 3(a) and (b) reveals the variability in the rate β1 and β2 at which individuals come into contact with each other, which affects the disease's spread.Fig. 3Random noise occurring in the parameters for formulation of Brownian Motion in the influenza disease model, where (a) noise in β1 (b) noise in β2 (c) noise in α (d) noise in δ (e) noise in γ and (e) noise in γ1, where diffusion rate σ = 0.1.Fig. 3
This approach takes into account the inherent uncertainty and randomness in real-world interactions. By incorporating noise and random fluctuations in the contact rate parameters, the model reflects the unpredictable nature of human behavior and social dynamics. It is noticed that, in between weeks 25 to 80 contact rate vary rapidly and in week 80 it strikes its maximum level. That indicates the high transmission of disease from exposed to susceptible. We also see that, transmission rate β2 fluctuates more within time interval 20 to 40 weeks. On week 30, it approaches to maximum level, that indicated high probability of outbreak for transmission of disease from infected to susceptible population.
Fig. 3(c) illustrates the variation in the rate at which the illness spreads from an infected person to vulnerable people. This approach acknowledges the uncertainty and randomness associated with the transmission of the disease. By incorporating noise and random fluctuations in the infection rate parameter, the model captures the unpredictability of how contagious the disease is in different situations. It considers factors such as variations in individual susceptibility, changes in behaviour, and external influences that can affect the likelihood of transmission. It is visible that, along week 20 to 40 the fluctuation og α lies in meridian level, but along 60 to 120 weeks the randomness becomes high. That means, on that time interval disease persists in the region.
By analyzing the fluctuation in death rate from Fig. 3(d), we observe that, by the progression of time, after 20 weeks, the death rate goes under control. This means, disease outbreak probability mitigated. Fig. 3(e) and (f) shows how these differences affect the model's forecast accuracy and dependability. The uncertainty or variability in the rate at which people recover from the illness is referred to as the noise in the parameter recovery rate. Similarly, the noise in treatment rate represents the uncertainty in the effectiveness or availability of treatments. It is observed that, recovery rate lies always in satisfactory level with the density of expose and infected individuals. Along weeks 60 to 100, the recovery rate become high, that indicates the seasonal variation when disease transmitting probability reduced rapidly. On the other hand, it is visualized that, the treatment effectiveness lies in medium level along weeks 0 to 40. After 60 weeks, the effectiveness increases as more population recovers by natural immunity or seasonal variation. The analysis, supports that γ1 can play significant role to mitigate the rates β1, β2 and α to reduce disease burden.
The qualitative behaviour of the stochastic model, which aids in our understanding of the persistence and extinction of a disease, is examined in this section. In stochastic models, randomness plays a significant role, and studying qualitative behaviour provides insights beyond average outcomes. We can assess whether a disease will persist in the population or eventually go extinct by looking at its extinction and persistence. This approach aids in evaluating the disease's long-term effects and the efficacy of control strategies. We can locate important thresholds or tipping points that decide whether a disease will persist or vanish using qualitative behaviour analysis. These thresholds might be influenced by factors including population size, transmission rate, and recovery rate. Additionally, investigating qualitative behaviour enables us to comprehend how random fluctuations affect the dynamics of disease. Qualitative analysis aids in the identification of the potential consequences of these stochastic fluctuations because random events can result in unanticipated outcomes (Rao et al., 2012; Yang, 2016; Zhang et al., 2020).
We begin by introducing the following (a)Define R+d={χi∈Rd:χi>0,1≤i≤d}.(b)Let Ω,F,{Ft}t≥0,P be a complete probability space equipped with a filtration {Ft}t≥0 satisfying the usual conditions.
Consider, in general, an n-dimensional stochastic differential equation (SDE) of the form(4.1)dy(t)=F(y(t),t)dt+G(y(t),t)dB(t),t≥0,with initial condition y(t0)=y0∈Rd.
We define the differential operator L associated with equation (4.1) asL=∂∂t+∑i=1dFi(y,t)∂∂yi+12∑i,j=1dGT(y,t)G(y,t)ij∂2∂yi∂yj.
If the operator L acts on a sufficiently smooth function V:Rd×R+→R+, thenLV(y,t)=∂V∂t(y,t)+∇yV(y,t)⋅F(y,t)+12traceGT(y,t)∇y2V(y,t)G(y,t),where ∇yV and ∇y2V denote the gradient and Hessian matrix of V with respect to y, respectively.Theorem 1(Elhiwi, 2021; Yang, 2016) Let (S(0),E(0),I(0),R(0),T(0))∈R+5 be the initial conditions. Then, the system (2.1) admits a unique positive solution (S(t), E(t), I(t), R(t), T(t)) for all t ≥ 0. Moreover, the solution remains in R+5 with probability one.
By deriving the expression for the threshold quantity R0 associated with the model (1), we conclude the (a)If 0<R0<1, the disease-free equilibrium Λμ,0,0,0,0 is globally asymptotically stable. However, it becomes unstable when R0>1.(b)If R0>1, the endemic equilibrium (S∗, E∗, I∗, R∗, T∗) of system (1) is globally asymptotically stable (Rao et al., 2012; Yang, 2016).
Recently, several researchers have incorporated parameter perturbations into epidemic models to explore their dynamic behaviors (Arruda et al., 2021; Rao et al., 2012; Wang et al., 2017).
In this work, accounting for the effects of a randomly fluctuating environment, we introduce white noise perturbations into each equation of the model (1). It is assumed that environmental variability predominantly influences the parameters γ1 and γ2, which are modified as γi⟶γi+σiB˙i(t),i=1,2,where Bi(t) (i = 1, 2) are independent standard Brownian motions with Bi(0) = 0, and σi2(i=1,2) represent the intensities of the respective white noise terms.
Accordingly, the stochastic version of the deterministic model (2.1) can be formulated as (4.2)dS=[Λ−(β1E+β2I)S−μS]dt+σ1SdB1(t)dE=[(β1E+β2I)S−(α+μ)E]dt+σ2EdB2(t)dI=[αE−(μ+δ+γ+γ1)I]dt+σ2IdB2(t)dR=[γI−μR]dt+σ1RdB1(t)dT=[γ1I−μT]dt+σ1TdB1(t)From (4.2), the equation for the overall population size is derived as dN≤(Λ−μN)dtIt follows that,(4.3)limt→∞(S(t),E(t),I(t),R(t),T(t))≤ΛμBy (4.3), we takeS(t)=Λμ−E(t)−I(t)−R(t)−T(t)and substitute it into the second equation of model (4.2), consequently we can obtain from the following (4.4)dE=(β1E+β2I)Λμ−E−I−R−T−(α+μ)Edt+σ2EdB2(t)dI=[αE−(μ+δ+γ+γ1)I]dt+σ2IdB2(t)Given an initial condition E(0),I(0)∈R+2, it can be easily verified through straightforward calculations that the model (3) admits a unique disease-free equilibrium, denoted byP0=Λμ,0,0,0,0.
Throughout this work, we assume thatΩ,F,{Ft}t∈R,Pis a complete probability space equipped with a filtration {Ft}t∈R satisfying the usual conditions (i.e., the filtration is right-continuous, non-decreasing, and F0 contains all P-null sets).
Define the domainΓ=(E,I)∈R+2:0<E(t)+I(t)<Λμ.
The following lemma is quoted from (Allen, 2010; Yang, 2016) where it was proved and applied. We proceed as the similar patterns in the problem.Lemma 1(Yang, 2016). Let, x ∈ C[Ω × [0, ∞], (0, ∞)]. If there exist positive constants λ, μ such thatlogx(t)≥λt−μ∫0tx(s)ds+F(t),asymptotically stablefor all t ≥ 0, where F∈C[Ω×[0,∞],R] and limt→∞F(t)t=0, thenlim inft→∞1t∫0tx(s)ds≥Λμ,asymptotically stable
In the study of the dynamical behavior of population models, a fundamental requirement is to verify that solutions remain positive and exist globally. Motivated by the approach in (Yang, 2016), we first establish the global existence of solutions for the model (4.4). Based on this, we present the following Lemma 2(Yang, 2016) Suppose that (E(0), I(0)) ∈ Γ. Then, the subsystem (4.4) admits a unique solution (E(t), I(t)) for all t ≥ 0, and the solution remains within the set Γ with probability one. That is, Γ is almost surely a positively invariant set for the model (4.4).
In this section, we present a theorem that establishes conditions under which the equilibrium point of the model (2.1) is almost surely exponentially stable. The formulation of this theorem is inspired by the work in (Yang, 2016).
Let us define σ≔ min{σ1, σ2} and X(t)≔(E(t), I(t)). Before proceeding, we first state a property concerning the disease-free dynamics (i.e., when I = 0), as discussed in (Rao et al., 2012; Wang et al., 2017; Yang, 2016).Theorem 2If
R0=S0αβ2+S0β1(μ+δ+γ+γ1)(α+μ)(μ+δ+γ+γ1)<1
then the disease-free equilibrium (0, 0) of (1) is almost surely exponentially stable, i.e. the disease will become extinct with a probability of one.ProofLet
a1
be any fixed positive real number. We define the following stochastic process:(4.5)z(X(t))=a1E(t)+I(t),where clearly z(X(t)) > 0 for all t > 0.
Using this, we introduce a C^2^-class function V:R+2→R+ defined byV(X(t))=lnz(X(t)).
Applying Itô’s formula, the stochastic process V(X(t)) evolves according toV(X(t))=V(X(0))+∫0tLV(X(τ))dτ+G(t),where G(t) represents the local martingale term associated with the stochastic integral.G(t)=∫0t−a1σ2E(τ)z(X(τ))dB1(τ)+∫0t−σ2I(τ)z(X(τ))dB2(τ)and,LV(X(t))=1za1(β1E+β2I)Λμ−E−I−R−T−a1(α+μ)E+αE−(μ+δ+γ+γ1)I−1z2a12σ22E22+σ22I22We have the following regarding the quadratic variations of the stochastic integral G(t):∫0t−a1σ2E(τ)z(X(τ))2dτ≤σ22t,∫0t−σ2I(τ)z(X(τ))2dτ≤σ22tAccording to the strong law of large number of martingles (Srivastav et al., 2022; Yang, 2016), we have,limt→∞G(t)t=0asymptotically stableNext, we aim to demonstrate that LV(X(t))<0, establishing the asymptotic stability of the system. To facilitate this, we introduce the following normalized v(t)=E(t)z(X(t)),w(t)=I(t)z(X(t)).It is clear that for all t > 0,0<v(t)<1a1,0<w(t)<1,implying that the stochastic processes v(t) and w(t) are non-negative and bounded above by max1a1,1.
From the definition of z(X(t)) in (6), it follows thata1v(t)+w(t)=1for all t>0.Using this relation, we can now express the differential operator LV(X(t)) as LV(X)=a1(β1v+β2w)Λμ−E−I−R−T−(α+μ)v+(αv−(μ+δ+γ+γ1)w)−a12σ22v22−σ22w22≤a1β1vΛμ−E−I−R−T+β2wΛμ−E−I−R−T−(α+μ)v+[αv−(μ+δ+γ+γ1)w]−σ222[a12v2+w2]≤a1β1Λμ−E−I−R−T−a1(α+μ)+a1αv+a1β2Λμ−E−I−R−T−(μ+δ+γ+γ1)w−σ22[a12v2+w2]In a view of a1v + w = 1, we have from the following expression that,−σ22(a1v2+w2)≤−σ24(a1v+w)2=−σ24(a1v+w)It follows that,(4.6)LV(X)≤A1v+A2wwhere,A1=a1β1Λμ−E−I−R−T−a1(α+μ)+a1α−a1σ24A2=a1β2Λμ−E−I−R−T−(μ+δ+γ+γ1)−σ24Let, a1=β1Λμβ1+(α+μ)+σ24, then A1 = 0. The condition R0<1 is comparable to the inequality that β1Λβ1+(α+μ)+σ24≤(α+μ)μ+μσ24−αμΛNow let a2=β2Λμβ2+(μ+δ+γ+γ1)+σ24, which implies that A2 < 0.
As a result, we have 0 < w < 1, and LV(X(t))<0. Finally, by dividing both sides of (4.6) by t and taking the limit as t → ∞, we obtainlimt→∞sup1tlnzX(t)=lim supt→∞1t∫0tLV(X(τ))dτ<0asymptotically stable.This completes the required assertion. □
The stochastic EE, or persistence of E and I, of (5), a subsystem of (4.2), will be examined in this section under specific parametric constraints (Rao et al., 2012; Wang et al., 2017; Yang, 2016). Theorem 3. If R0=S0αβ2+S0β1(μ+δ+γ+γ1)(α+μ)(μ+δ+γ+γ1)>1 and R0∗=α(β1+β2)Λβ1+α+μ+σ222μ(μ+δ+γ+γ1)+μσ222−αΛ>1 for any initial value (E(0), I(0)) ∈ Γ of model (5). has the following lim inft→∞1t∫0tE(s)ds≥αβ1Λμβ1+α+μ+σ222μ(μ+δ+γ+γ1)+μσ222−αβ2Λ(R0∗−1)and,lim inft→∞1t∫0tI(s)ds≥μ(μ+δ+γ+γ1)+μσ222−αβ2Λ(R0∗−1)This suggests that the model (5) solutions show a high degree of persistence in the mean.ProofAn integration of the primary equation (perturbed equation) of (5) of the model (3) provides,β1Λμ⟨E(t)⟩−β1⟨E(t)2⟩−β1⟨E(t)I(t)⟩−β1⟨E(t)R(t)⟩−β1⟨E(t)T(t)⟩+β2Λμ⟨I(t)⟩−β2⟨E(t)I(t)⟩−β2⟨I(t)2⟩−β2⟨I(t)R(t)⟩−β2⟨I(t)T(t)⟩−(α+μ)⟨E(t)⟩=E(t)−E(0)t+σ2t∫0tE(s)dB1(s)We compute that,(4.7)⟨E(t)⟩=β1αΛ2μ2(β1+α+μ)⟨E(t)⟩−β1αΛμ(β1+α+μ)⟨E(t)2⟩(4.8)−β1αΛμ(β1+α+μ)⟨E(t)I(t)+E(t)R(t)+E(t)T(t)⟩+β1β2Λ2αμ2(β1+α+μ)⟨I(t)⟩−β1β2Λαμ(β1+α+μ)⟨E(t)I(t)+I(t)2+I(t)R(t)+I(t)T(t)⟩−β1α(α+μ)Λμ(β1+α+μ)⟨E(t)⟩−φ(t)≤β1β2Λαμ(β1+α+μ)⟨I(T)2+I(t)R(t)+I(t)T(t)⟩+β1β2Λ2αμ2(β1+α+μ)⟨I(t)⟩−φ(t)where, φ(t)=1β1+α+μE(t)−E(0)t+σ2t∫0tE(s)dB1(s)
Since, E(t),I(t)<Λμ, by using the strong law of large numbers for martingles (Arruda et al., 2021; Yang, 2016), we havelimt→∞E(t)t=0,limt→∞1t∫0tE(s)dB1(s)=0asymptotically stableObviously, lim~t→∞φ(t) = 0 is asymptotically stable. By applying Itô’s formula to the transformed equation (5) of the model (3), we (4.9)dlnμΛE≤β2αΛμIE−β2αΛμI−β2αΛμI2E−β2αΛμIRE−β2αΛμITE−(α+μ)+12σ22+σ2dB2(t)An integration of (10) αβ1Λμβ1+α+μ+σ222IE=1+1μβ1+α+μ+σ222αβ1I2E+⟨I⟩+σ2dB2(t)t+lnμΛ(E(t)−E(0))tAs long as since 0<E(t)<Λμ, it follows that −∞<lnμΛE(t)<0. For any 0 < ϵ < 1, there exist T = T(w) > 0 and a set Ωϵ~ such that P(Ω~ϵ) ≥ 1 − ϵ. For all t ≥ T(w) and w ∈ Ωϵ~,αβ1Λμβ1+α+μ+σ222IE<1,such that,αβ1Λμβ1+α+μ+σ222∫0tI(s)E(s)ds<tSince 0<E(t),I(t)<Λμ, it follows (4.10)αβ1Λμβ1+α+μ+σ222I(t)<E(t)Thus, by applying Itô’s formula to the transformed equation (5) of the model (3), we (4.11)d(lnI)=αβ2Λμ+αEI−(μ+δ+γ+γ1)+12σ22+σ2dB2(t)Integrating both sides from 0 to t, we lnI(t)t≥αβ2Λμβ2+μ+δ+γ+γ1+σ222+αβ2Λμ−(μ+δ+γ+γ1)+12σ22+α2β22Λμ2β2+μ+δ+γ+γ1+σ222+αβ2[I(t)]+φ(t)+σ2dB2(t)t+lnI(0)tSince,limt→∞φ(t)+σ2dB2(t)t+lnI(0)t=0by Lemma 1, it follows limt→∞inf(I(t))≥αβ2Λμβ2+μ+δ+γ+γ1+σ222+αβ2Λμ−(μ+δ+γ+γ1)+12σ22≥μ(μ+δ+γ+γ1)+μσ222−αβ2Λ(R0−1)>0Finally, from the last equality in (11), we limt→∞infE(t)>αβ1Λμβ1+α+μ+σ222μ(μ+δ+γ+γ1)+μσ222−αβ2Λ(R0∗−1)also this yeilds,limt→∞inf(I(t))>0asymptotically stableThis completes the required proof.
In the following sections, we present several theorems regarding the extinction and persistence of influenza disease in a stochastic sense. To this end, a mathematical model for influenza with random perturbations is formulated as follows (Gray et al., 2011; Rao et al., 2012):(4.12)dSdt=Λ−β1E(t)S(t)−β2S(t)I(t)−μS(t)−ρS(t)I(t)dB(t)−ρS(t)E(t)dB(t)dEdt=β1E(t)S(t)+β2S(t)I(t)−(α+μ)E(t)+ρS(t)E(t)dB(t)dIdt=αE(t)−(μ+δ+γ+γ1)I(t)+ρS(t)I(t)dB(t)−ρ(R(t)+T(t))dB(t)dRdt=γI(t)−μR(t)+ρR(t)dB(t)dTdt=γ1I(t)−μT(t)+ρT(t)dB(t)where B(t) represents the standard Brownian motion, with ρ^2^ > 0 and the intensity of the white noise.
In this section, we explore the conditions for the extinction of influenza spread, as studied in (Rao et al., 2012; Srivastav et al., 2022; Wang et al., 2017; Zhang et al., 2020). Here, we establish the following ⟨y(t)⟩=1t∫0ty(s)dsand we have calculated the threshold quantity. R0=S0αβ2+S0β1(μ+δ+γ+γ1)(α+μ)(μ+δ+γ+γ1)
and let, R~0=(β1+β2)Λμ1(α+μ)+12ρ2Λμ
A useful lemma relevant to this study is stated Lemma 3(Bayram et al., 2018; Wang et al., 2017) Let
M={Mt}t≥0
be a real-valued, continuous local martingale satisfying
M0 = 0.Then, the following •If lim~t→∞⟨M,M⟩t~ = ∞, it follows thatlimt→∞Mt⟨M,M⟩t=0.•Moreover, iflimt→∞sup⟨M,M⟩tt<∞,thenlimt→∞Mtt=0.Theorem 4Let, the solution of system (13) be (S(t), E(t), I(t), R(t), T(t)), with initial value
(S(0),E(0),I(0),R(0),T(0))∈R+5. If(1)ρ2>max(β1+β2)22(μ+δ+γ+γ1),(β1+β2)μΛ or,(2)R0<1, R~0<1 and ρ2≤(β1+β2)μΛ
Then,(4.13)lim supt→∞logI(t)t≤−(μ+δ+γ+γ1)+(β1+β2)2ρ2<0is asymptotically stable if condition (1) holds.(4.14)lim supt→∞logI(t)t≤(β1+β2)Λμ1−1R~0<0is asymptotically stable if condition (2) holds. In addition,limt→∞S(t)=Λμ=S0,limt→∞E(t)=0,limt→∞I(t)=0,limt→∞R(t)=0,limt→∞T(t)=0is asymptotically stable. In other words, the disease I is stochastically extinct exponentially with probability one.ProofBy integrating system (13), we S(t)−S(0)t=Λ−β1⟨S(t)E(t)⟩−β2⟨S(t)I(t)⟩−μ⟨S(t)⟩−ρ⟨S(t)I(t)⟩dB(t)E(t)−E(0)t=β1⟨S(t)E(t)⟩+β2⟨S(t)I(t)⟩−(α+μ)⟨E(t)⟩+ρ⟨S(t)E(t)⟩dB(t)I(t)−I(0)t=α⟨E(t)⟩−(μ+δ+γ+γ1)⟨I(t)⟩+ρ⟨S(t)I(t)⟩dB(t)−ρ⟨R(t)+T(t)⟩dB(t)R(t)−R(0)t=γ⟨I(t)⟩−μ⟨R(t)⟩+ρ⟨R(t)⟩T(t)−T(0)t=γ1⟨I(t)⟩−μ⟨R(t)⟩+ρ⟨T(t)⟩dB(t)Then we have,S(t)−S(0)t+E(t)−E(0)t+I(t)−I(0)t+R(t)−R(0)t+T(t)−T(0)t≤Λ−μ⟨S(t)⟩−(μ+δ+γ+γ1)−αγμ+α⟨I(t)⟩≤Λ−μ⟨S(t)⟩−(μ+α)(μ+δ+γ+γ1)−αγ(μ+α)⟨I(t)⟩Therefore,⟨S(t)⟩=−1μS(t)−S(0)t+E(t)−E(0)t+I(t)−I(0)t+R(t)−R(0)t+T(t)−T(0)t+Λμ−1μ(μ+α)(μ+δ+γ+γ1)−αγ(μ+α)⟨I(t)⟩By applying(4.15)limt→∞⟨S(t)⟩=Λμ−1μ(μ+α)(μ+δ+γ+γ1)−αγ(μ+α)⟨I(t)⟩(4.16)dlogI(t)≤(β1+β2)S−(μ+δ+γ+γ1)−12ρ2S2dt+ρSdB(t)(4.17)logI(t)−logI(0)t≤(β1+β2)⟨S(t)⟩−(μ+δ+γ+γ1)−12ρ2⟨S(t)2⟩+12∫0tS(r)dB(r)Let, β = (β1 + β2)
By putting the value of S(t) from equation (16) we have,logI(t)−logI(0)t≤βΛμ−1μ(μ+α)(μ+δ+γ+γ1)−αγ(μ+α)⟨I(t)⟩−(μ+δ+γ+γ1)−12ρ2Λμ−1μ(μ+α)(μ+δ+γ+γ1)−αγ(α+μ)⟨I(t)⟩2+ρt∫0tS(r)dB(r)≤βΛμ−μ+δ+γ+γ1μ+α⟨I(t)⟩−(μ+δ+γ+γ1)−12ρ2Λμ2−μ+δ+γ+γ1μ+α2⟨I(t)⟩2+2Λμμ+δ+γ+γ1μ+α⟨I(t)⟩+ρt∫0tS(r)dB(r)From this we have,logI(t)−logI(0)t≤βΛμ−(μ+δ+γ+γ1)−12ρ2Λμ2−β(μ+δ+γ+γ1)(α+μ)⟨I(t)⟩+2Λμμ+δ+γ+γ1α+μ⟨I(t)⟩−12ρ2−μ+δ+γ+γ1α+μ2⟨I(t)⟩2+ρt∫0tS(r)dB(r)≤βΛμ−(μ+δ+γ+γ1)+12ρ2Λμ2−β(μ+δ+γ+γ1)α+μ⟨I(t)⟩+2Λμμ+δ+γ+γ1α+μ⟨I(t)⟩−12ρ2−μ+δ+γ+γ1α+μ2⟨I(t)⟩2+ρt∫0tS(r)dB(r)Consequently, further calculation gives,logI(t)−logI(0)t≤βΛμ1−μ(μ+δ+γ+γ1)+12ρ2Λμ2βΛ−β(μ+δ+γ+γ1)μ+α⟨I(t)⟩+2Λμμ+δ+γ+γ1μ+α⟨I(t)⟩−12ρ2−μ+δ+γ+γ1μ+α2⟨I(t)⟩2+ρt∫0tS(r)dB(r)≤βΛμ1−1R~−β(μ+δ+γ+γ1)(μ+α)⟨I(t)⟩+2Λμμ+δ+γ+γ1μ+α⟨I(t)⟩−12ρ2−μ+δ+γ+γ1μ+α2⟨I(t)⟩2+ρt∫0tS(r)dB(r)If condition (2) holds true, thenlim supt→∞logI(t)t≤(β1+β2)Λμ1−1R~<0and conclusion (15) is proved. Next, according to inequality (18) we have,logI(t)−logI(0)t≤(β1+β2)⟨S(t)⟩−(μ+δ+γ+γ1)−12ρ2⟨S(t)⟩2+ρt∫0tS(r)dB(r)≤−12ρ2S(t)−βρ2+β2ρ2−(μ+δ+γ+γ1)+ρt∫0tS(r)dB(r)If condition (1) holds true, thenlogI(t)t≤β2ρ2−(μ+δ+γ+γ1)+ρt∫0tS(r)dB(r)+logI(0)tand conclusion (15) is proved. We have,limt→∞logI(t)t≤−(μ+δ+γ+γ1)+β2ρ2<0is asymptotically stable.According to (14) and (15) we have,limt→∞I(t)=0Now, from the equations of system (13), it follows that,R(t)=e−μtR(0)+∫0tγI(r)eμrdtBy applying the L'Hospitals rule to the previous result, we havelimt→∞R(t)=0Now, adding all the equation of system (13), it follows that,N(t)≤e−μtN(0)+Λμe−μt⇒S(t)+E(t)+I(t)+R(t)+T(t)≤S(0)+E(0)+I(0)+R(0)+T(0)+Λμe−μteμt⇒limt→∞S(t)=limt→∞S(0)+E(0)+I(0)+R(0)+T(0)+Λμe−μteμt−(E(t)+I(t)+R(t)+T(t))It follows that,limt→∞S(t)=ΛμTherefore, if condition (2) holds and for R0<1 and R~<1, the disease induced density I will be gradually extinct from the community with probability one. Thus, we have concluded the proof.
This section concerns the persistence of disease of system (13).
Theorem 5. Suppose that μ>ρ22. Let, (S(t), E(t), I(t), R(t), T(t)) be any solution of model (13) with initial conditions (S(0),E(0),I(0),R(0),T(0))∈R+5. If R0>1 and R~>1, thenlimt→∞⟨S(t)⟩=Λμ−1μ(μ+α)(μ+δ+γ+γ1)−αγμ+α⟨I(t)⟩=Λμ−βΛμ1−1Rβ−2Λμlimt→∞⟨I(t)⟩=βΛμ1−1Rμ+δ+γ+γ1α+μβ−2Λμlimt→∞⟨R(t)⟩=γ(μ+δ+γ+γ1)βΛμ1−1R~β−2ΛμThis implies that the disease I in the model (4.12) is stochastically persistent on average.
In this section, we analyze the persistence and extinction of disease dynamics in the exposed and infected compartments. For this purpose we have assumed the initial conditions same as SDE model analysis, and the contact rates β1 = 0.35, β2 = 0.35. The diffusion rate to each compartment is considered σ = 0.1. Fig. 4(a) illustrates the scenario in the exposed compartment by varying the contact rates β1 and β2, which represent the average rate at which susceptible individuals interact with exposed and infectious individuals, facilitating potential exposure and disease transmission (Srivastav et al., 2022; Zhang et al., 2020). We observe that when the contact rates are high (β1 = 0.45, β2 = 0.45), the disease is more likely to spread quickly through the population, resulting in a larger number of individuals becoming exposed within a shorter time frame. Consequently, within 0 to 20 weeks, exposed density strikes at 360. After that it gradually decreases but fluctuates between 50 and 100. That indicates, disease will persist for a long time. With a moderate contact rate, the disease may still spread through the population, but at a slower rate compared to the high contact rate scenario. When the contact rate is low (β1 = 0.25, β2 = 0.25), disease transmission is limited, and the spread of the disease occurs at a slower pace. This can result in a lower number of individuals being exposed, leading to a decrease in the persistence of the disease in the exposed compartment. In this case, exposed density strikes to peak level 180 in week 23. After that the convergence to the zero level visualized. This fluctuation always lower than the previous one.Fig. 4Persistence of influenza disease model with random noise, where (a) influence of β1 and β2 in exposed (b) influence of β1 and β2 in infected (c) influence of γ in infected (d) influence of γ1 in infected, where, β1 = 0.35, β2 = 0.35, α = 0.25, γ = 0.08, γ1 = 0.05, μ = 0.01, δ = 0.01, and diffusion rate σ = 0.1.Fig. 4
Fig. 4(b) expresses that, contact rates β1 and β2 has positive influence on infected population density. For the case when β1 = 0.25, β2 = 0.25 infected individuals approaches to maximum level within weeks 37. On the other hand, when β1 = 0.45, β2 = 0.45, then a rapid increase in infected compartment noticed. Within 22 weeks, infection strikes peak level 510. Analyzing the levels of these two curves, we realized that the higher the contact rates, the outbreak of disease spreading from infected compartment rises. Thus red curve indicates, disease will persist for longer period.
Fig. 4(c) shows that the recovery rate (γ) in the infected compartment plays a crucial role in determining whether a disease persists or goes extinct. We observe that, a high recovery rate (γ = 0.08) represent that infected individuals have a shorter duration of infection before they recover and become immune or no longer infectious. In this scenario, the disease is more likely to be contained and eventually extinguished. Within 0 to 25 weeks, the density of infection gradually increase and at week 25 the maximum infection case reported. As infected individuals recover quickly, they are removed from the infected compartment, reducing the trend of infectious individuals and limiting further transmission. Moreover, after 40 weeks the random movement in infection density fluctuates around level 100 to 110. That indicates, disease will extinct within 60 weeks. On the other hand, low recovery rate (γ = 0.02) implies that infected individuals take a longer time to recover or have a higher probability of severe illness or death. In this case, the disease is more likely to persist and potentially spread within the population. The longer duration of infection allows for a greater opportunity for transmission to susceptible individuals, leading to a higher likelihood of sustained transmission and a slower decline in the infected compartment. For this case, it is visualized that, the peak level of infection increases up to 500 in week 25. After 25 weeks, the infection rate randomly moves around 200 to 250. That indicates, disease will persist for long time.
Fig. 4(d) shows that, the treatment rate (γ1) can have a significant impact in the infected compartment on the persistence or extinction of a disease. If the treatment rate absent i.e., (γ1 = 0), the infected individuals may not receive effective treatment in a timely manner. This can lead to the continued transmission of the disease within the population, resulting in its persistence. In this scenario, infected individuals may not fully recover or may experience relapses, allowing the disease to persist and potentially spread to new individuals. We observe that, within 20 to 40 weeks the rapid increase in infection class occurs. After, 80 weeks the outbreak still rises in community. Conversely, if the treatment rates are high (γ1 = 0.05) then the infected individuals receive prompt and effective treatment, the disease can be brought under control and eventually eliminated. Adequate treatment can reduce the infectiousness of individuals, shorten the duration of their illness, and prevent transmission to others. With a high treatment rate, the number of infected individuals decreases over time, and if this trend continues, the disease may eventually be eradicated from the population. It is noticed that, after striking to peak level, the infection lessens between weeks 40 to 60 For this case, the red curve always stays below the blue curve, which ensures the reduction of infection density satisfactorily.
The capacity of viruses to survive, host behavior, and immune system performance are all impacted by seasonal fluctuations in temperature, humidity, and precipitation. These factors ultimately affect the effectiveness of transmission. The geographical distribution of sickness in human populations and in populations of domestic and wild animals is also influenced by the seasonal migration of susceptible and afflicted hosts. Numerous studies have explored the combined effects of viral infections enhanced by cytokines and cell-free infections, proposing a partial differential equation (PDE) model that incorporates geographic variability. In these studies, the basic reproduction number R0 was defined as the spectral radius of the sum of two linear operators associated with viral infections, one of which is cytokine-enhanced and the other is cell-free. These studies also confirmed the threshold behavior of this number. Recent works have examined an efficient optimized hybrid block approach designed to integrate generic second-order initial value problems (IVPs) of ordinary differential equations (ODEs). By introducing an adaptive step-size formulation, researchers proposed an enhanced method, which they found to be a suitable alternative to existing solvers, offering similar features. Furthermore, a SIR (Susceptible-Infected-Recovered) epidemic model with free boundaries, non-local infection, and diffusion was proposed, and the existence of a solitary global solution was established. Additionally, we have developed and studied an SIR model where incidence rates are represented by a Monod-type equation and the recovery process is non-linear. Recent research has also introduced a non-standard finite difference (NSFD) scheme for the model, with the denominator function selected to ensure the boundedness of the solutions in the proposed scheme. Researchers who have studied the sign of wave speed for disable traveling waves to a two-species competitive system have lately scientifically simulated the dynamics of two species vying for a shared scheme (Gray et al., 2011; Nipa, 2020).
Seasonal contact and other seasonal behaviors can be captured by stochastic or deterministic epidemic models using periodic coefficients. Predicting the progression of disease outbreaks in periodic environments presents greater complexity than in continuous conditions. Recent research has focused on epidemic models in both deterministic and stochastic periodic settings. To better understand the influence of seasonality and demographic variation on the probability of disease outbreaks, we extend some of these stochastic studies to incorporate a non-homogeneous stochastic model. This model specifically addresses the transmission of diseases in periodic environments, taking into account the variations in both time and population dynamics that characterize the seasonal fluctuations and demographic changes that affect disease spread (Arruda et al., 2021; Niu et al., 2021).
In specifically, the probability of an epidemic of a disease is quantified using an approximate multi-type branching process. A system of ordinary differential equations derived from the inverted Kolmogorov differential equations provides an approximation for the likelihood of a disease outbreak. The analysis reveals that the probability of an outbreak follows a cyclical pattern, which is strongly influenced by the precise timing at which an infected or exposed individual is introduced into a population of entirely susceptible individuals or susceptible vectors. Numerical simulations, using periodic transmission rates for various compartments, show that the peak transmission rates from the exposed and infected compartments to the susceptible compartments coincide with the moments when the probability of an epidemic is at its highest.The main objective of this section 1.Observing how the infectious illness transmission model's dynamics are affected by seasonal fluctuation.2.Introducing and briefly discussing the approximation of the branching process.3.Discovered the likelihood of an epidemic in relation to a seasonal and deterministic environment.The findings of this section regarding the goal 1.Seasonal fluctuation affects the number of populations exposed to and afflicted with infectious illnesses.2.The deterministic and periodic settings have different basic reproduction numbers.3.The likelihood of a disease outbreak is significantly affected by the initial number of hosts in the population.4.The probability of disease outbreaks is significantly affected by seasonal variations.
In this section, we model the transmission of influenza using a system of non-linear ordinary differential equations. The total human population, denoted as N, is divided into three susceptible individuals, represented by S(t), exposed individuals, denoted by E(t), and infectious individuals, represented by I(t). That means the first three equations of (1) has been considered. We assume that there exist disease-related deaths. the susceptible humans get the infection through contact with the exposed and infectious humans. The human population can be resented as SEIS model, with the total population density expressed as N = S + E + I. The human demography is modeled using the parameter Λ, which represents the recruitment rate of healthy hosts. The parameters μ and δ denote the natural death rate for all compartments and the disease-induced death rate, respectively. Additionally, γ and γ1 represent the recovery rate from infection and the rate of successful progression to the treated compartment, respectively. After a latent period, the exposed population becomes infected at a rate α. Furthermore, the transmission rate from exposed individuals to susceptible individuals is given by β1(t), while the transmission rate from infected individuals to susceptible individuals is represented by β2(t). Both β1(t) and β2(t) are periodic functions of time with the same period. That means,(5.1)βi(t)=βi(t+m),i=1,2,t∈(−∞,∞)with period m > 0. Here, i = 1, 2 denotes the human compartments. The host (human) population size is constant, N(t) and therefore at disease free equilibrium state, S(t)=N(t)=Λμ and rest all E = I = R = T = 0. This the disease free equilibrium (DFE) for human is, D=(S0,E0,I0,R0,R0)=Λμ,0,0,0,0. The existence and stability of a periodic endemic equilibrium could also be examined according (Srivastav et al., 2022; Wang et al., 2017; Zhang et al., 2020). In particular, we focus on the dynamics near the disease-free equilibrium (DFE) to analyze the probability of an outbreak within a more general stochastic framework.
The following assumptions were introduced by Wang and Zhao (Allen, 2010; Nipa, 2020). The system of ordinary differential equations (ODE) can be expressed as the following dxidt=Fi(t,x)−Vi(t,x)≡fi(t,x),i=1,2,3,4,5where x = (x1, x2, x3) = (S, E, I) and Vi=Vi+−Vi−. The following assumptions are (M1)The mappings Fi(t,x), Vi+(t,x), and Vi−(t,x) are continuous and non-negative on R×R5+ and are C^1^ with respect to x.(M2)The mappings Fi(t,x), Vi+(t,x), and Vi−(t,x) are p-periodic in t, where p > 0.(M3)If xi = 0, then Vi−(t,x)=0.(M4)Fi(t,x)=0 for i = 1, 4, 5.(M5)If E = I = 0, then Fi(t,x)=Vi+(t,x)=0 for i = 2, 3.(M6)If E = I = 0, the disease-free equilibrium (DFE) is stable.(M7)The spectral radius of the matrix FV−1 is less than one.
In the special case where the transmission rates are constant, β1=β1~ and β2=β2~, it is straightforward to determine the basic reproduction number using the next-generation matrix, as described in Section A. This results R0~=β1S0(μ+δ+γ+γ1)+β2S0α(α+μ)(μ+δ+γ+γ1).An epidemic of illness happens when this threshold parameter is greater than 1. Type or target reproduction numbers are other variations of the basic reproduction number that have also been established. Furthermore, reducing the threshold quantity through control methods also lowers the likelihood of an epidemic. When modeling epidemics with time-periodic coefficients using non-autonomous systems of differential equations, developed sufficient conditions for existence of R0 (conditions (M1)-(M7)) (Allen, 2010; Nipa, 2020; Srivastav et al., 2022; Zhang et al., 2020). These conditions depend on the monodromy matrix for E and I (5.2)dEdt=(β1E+β2I)S−(α+μ)EdIdt=αE−(μ+δ+γ+γ1)IFor X=(E,I)T, the linear system can be represented as dXdt=F(t)−V(t)X, where,(5.3)F(t)=β1S0β2S000andV(t)=α+μ0−αμ+δ+γ+γ1The systems (1) and (20) satisfy the assumptions (M1)-(M7). Only in special cases, does there exist an explicit solution for R0~ and R0 (Allen, 2010). However, the basic reproduction number, R0, for the periodic system can be numerically estimated by analyzing the following system of differential dXdt=F(t)⋌−V(t)X.Specifically, if ΩF⋌−V(t) represents the monodromy matrix and the spectral radius satisfies the condition ρ(ΩF⋌−V(p))=1, then the value of ⋌ corresponds to R0 (Nipa, 2020).
To obtain a numerical approximation of R0, we employ an iterative approach by constructing an increasing sequence ⋌i, where i = 1, 2, …, ensuring that the difference between consecutive values satisfies the constraint ⋌i+1 − ⋌i < 0.001. This process continues until we reach a specific value ⋌i that leads to a non-trivial periodic solution. The periodic solution satisfies the condition |X(t)−X(t+p)|<ϵ for an arbitrarily small ϵ and sufficiently large t.
To determine R0, we implement this numerical approximation method, which was originally proposed by Posny and Wang (Nipa, 2020).
In the context of infectious disease transmission, we construct a Continuous-Time Markov Chain (CTMC) model, denoted as (2.1).
In this model, we assume that the state variables representing the number of individuals in different compartments are discrete-valued, while time is considered to be continuous, with t∈0,∞. Specifically, we define the variables as S(t),R(t),T(t)∈{0,1,2,3,…,N},andE(t),I(t)∈{0,1,2,3,…,N}.This means that the susceptible (S), recovered (R), and treated (T) compartments take discrete values ranging from 0 to N, where N represents the total population size. Similarly, the exposed (E) and infectious (I) compartments are also represented as discrete values within the same range.
The infinitesimal transition probabilities, which describe the likelihood of transitioning from one state to another over an infinitesimally small time interval, are systematically presented in Table 2. These transition probabilities form the foundation for modeling the stochastic dynamics of the infectious disease spread over time.Table 2Quantifying the probability of transitions between various states (such as susceptible, infected, or recovered) within a very short time span is the task of representing infinitesimal transition probabilities for the stochastic infectious transmission model.Table 2Event,iDescription****(ΔS, ΔE, ΔI, ΔR, ΔT)****Probabilities1Rate of healthy human birth(1, 0, 0, 0, 0)ΛΔt + O(Δt)2Rate of natural death of healthy human(−1, 0, 0, 0, 0)μS(t)Δt + O(Δt)3Rate of disease transmission from exposed to susceptible(−1, 1, 0, 0, 0)β1E(t)S(t)Δt + O(Δt)4Rate of disease transmission from infected to susceptible(−1, 0, 1, 0, 0)β2I(t)S(t)Δt + O(Δt)5Rate of disease progression from exposed to infected(0, − 1, 1, 0, 0)αE(t)Δt + O(Δt)6Rate of natural death of exposed human(0, − 1, 0, 0, 0)μE(t)Δt + O(Δt)7Rate of Disease induced death from infected human(0, 0, − 1, 0, 0)δI(t)Δt + O(Δt)8Rate of natural death of infected human(0, 0, − 1, 0, 0)μI(t)Δt + O(Δt)9Rate of recovery from infected human(0, 0, − 1, 1, 0)γR(t)Δt + O(Δt)10Rate of treatment from infected human(0, 0, − 1, 0, 1)γ1I(t)Δt + O(Δt)11Rate of natural death of recovered human(0, 0, 0, − 1, 0)μI(t)Δt + O(Δt)12Rate of natural death of treated human(0, 0, 0, 0 − 1)μT(t)Δt + O(Δt)13No change(0, 0, 0, 0, 0)1 − ∑(t)Δt + O(Δt)
We define the variations in the susceptible, exposed, infected, and treated populations over a small time interval as ΔS(t),ΔE(t),ΔI(t),ΔR(t),ΔT(t). For instance, the change in the infected population is given ΔI(t)=I(t+Δt)−I(t).Let ∑(t) represent the total transition rate, expressed ∑(t)=Λ+μS(t)+β1E(t)S(t)+β2I(t)S(t)+αE(t)+μE(t)+(μ+δ+γ+γ1)I(t)+μR(t)+μT(t).
Incorporating seasonal variations into the transmission rates introduces time dependence, making the process non-homogeneous.
Let the state vector be defined as W(t)=(S(t),E(t),I(t),R(t),T(t)), and let the change in state over a small time interval Δt be represented as ΔW(t)=W(t+Δt)−W(t). For example, the probability of disease transmission from the exposed compartment to the susceptible compartment during a short time interval Δt is given PΔW(t)=(−1,1,0,0,0)∣W(t)=β1S(t)E(t)Δt+O(Δt).
We use branching process theory and probability generating functions (pgfs) to approximate the likelihood of avoiding a large outbreak (minor outbreak) for the Continuous Time Markov Chain (CTMC) model near the disease-free equilibrium. When E(t) and I(t) reach zero, they remain in that state, marking an absorbing state where disease transmission ceases. Both stochastic and deterministic models predict that, as t → ∞, the number of infectious individuals tends to zero with E(t) → 0 and I(t) → 0. However, in the stochastic model, extinction occurs at finite times. We focus on the dynamics at the onset of an epidemic when most of the population is still susceptible (Allen, 2010; Nipa, 2020).
The branching process approximates the disease-free equilibrium and is linearized by the SIS and SEIS models. For a small number of initially infected individuals, the process either reaches zero or expands exponentially. The CTMC model's branching process approximation accounts for these two outcomes. A significant epidemic occurs when the number of infectious individuals grows, while a minor outbreak happens with only a small increase in new cases. If the susceptible population is sufficiently large, the branching process provides a good approximation of the CTMC model, distinguishing between a large epidemic and a minor outbreak.
The branching process models E(t) and I(t), ignoring S(t), R(t), T(t). Starting with a few infectious individuals, the branching process approximation follows the CTMC near the disease-free equilibrium, where rates are linear (see Table 2).
The branching process approximation is based on three (1)The behavior of each infectious individual is independent.(2)Each infectious individual has the same recovery and transmission probabilities.(3)The susceptible population is sufficiently large.
We use two probability generating functions (pgfs) to study extinction one for each infectious individual (offspring pgf) and one for the entire infectious class E(t) and I(t). The offspring pgf is most critical.
Finally, we approximate the non-homogeneous stochastic process near the Disease-Free Equilibrium (DFE), leading to a multi-type branching process for E(t) and I(t) with periodic parameters β1(t) and β2(t), defined by the events 3-12 in Table 2. Then, the transition probability from state (e, i) at time t to state (j, k) at time t + Δt is defined as (5.4)p(e,i),(j,k)(t,t+Δt)=P{(E(t+Δt),I(t+Δt))=(j,k)|(E(t),I(t))=(e,i)}or, it can be written p(e+j,i+k)(e,i)(Δt)=Prob{E(t+Δt),I(t+Δt)=(j,k)|(E(t),I(t))=(e,i)}again we can express this as p(s,i),(l,k)(t+Δt)=P{S(t+Δt),I(t+Δt)=(l,k)|(S(t),I(t))=(s,i)}where, t < t + Δt.
DefinedP1dt=−(β1(t)S0+μ+α)[fE(P1,P2,t)−P1]=F1(P1,P2,t),dP2dt=−(β2(t)+μ+δ+γ+γ1)[fI(P1,P2,t)−P2]=F2(P1,P2,t).These equations are solved backward in time. Let tk = −kΔt, where Δt > 0. For instance, Euler's method provides the Pj(tk+1)≈Pj(tk)−ΔtFj(P1(tk),P2(tk),tk),Pj(0)=0,j=1,2Therefore, for (E(t), I(t)) = (e, i), the probability of an outbreak is approximated (5.5)Poutbreak(t)=1−[P1((t+p)−t)]e[P2((t+p)−t)]i.
When the transmission rates are assumed to be constant—meaning the periodic transmission rates βj(t), j = 1, 2 are replaced by their respective average values—an explicit expression for the probability of an outbreak can be derived (Nipa, 2020; Zhang et al., 2020). This probability is formulated in terms of the average threshold factor R0, as (5.6)Poutbreak=1−[qE]e[qI]iwhere, the value of qE and qI are(5.7)qE=β1S0β1S0+μ+α1R02+μ+αβ1S0+μ+αqI=β2S0β2S0+μ+δ+γ+γ11R~02+μ+δ+γ+γ1β2S0+μ+δ+γ+γ1In this particular case, the stochastic process is modeled as a time-homogeneous Markov chain, meaning that the transition probabilities between states remain constant over time (Arruda et al., 2021; Gray et al., 2011; Nipa, 2020; Niu et al., 2021).
The results from the previous sections are demonstrated through several numerical examples. Moreover, in this study, five sets of periodic transmission rates are considered for parameters β1, β2 and γ. The other parameters average values are estimated μ=0.01,α=0.25,γ1=0.05,δ=0.01The periodic nature in transmission rates are (5.8)β1(t)=β1ˇ1+σ1cos2πtp,β2(t)=β2ˇ1+σ2cos2πtpβ1(t)=β1ˇ1+σ1cos2πtp,β2(t)=β2ˇ1+σ2sin2πtpβ1(t)=β1ˇ1+σ1sin2πtp,β2(t)=β2ˇ1+σ2cos2πtpβ1(t)=β1ˇ1+σ1sin2πtp,β2(t)=β2ˇ1+σ2sin2πtpβ1(t)=β1ˇ1+amp×sin2π(t−phi)period,β2(t)=β2ˇ1+amp×sin2π(t−phi)period
In the numerical examples, the following parameters are period p = 0.1, amplitudes σ1 = σ2 = 0.2, phase ϕ = 0.25, and average transmission rates β1ˇ=0.35 and β2ˇ=0.4. The simulation time frame is one year. Both natural and disease-related mortality rates are considered in the outbreak probability. A summary of the parameter values is provided in Table 3. These values, though hypothetical, are biologically relevant and serve as illustrative estimates. For example, a recovery rate γ = 0.1 implies an average recovery time of 5 units, or approximately 3 months.Table 3Detailed explanation of the model parameters and their descriptions.Table 3NotationDefinition**Value**SouceS(t)Number of total susceptible human[500, 550][5, 8]E(t)Number of total exposed human(Khanh, 2016; Rosyada & Hariyanto, 2019)[9, 23, 37]I(t)Number of total infected human(Khanh, 2016; Rosyada & Hariyanto, 2019)[12, 15]R(t)Number of total recovered human(Centers for Disease Control and; Rosyada & Hariyanto, 2019)[11, 16]T(t)Number of total treated human(Centers for Disease Control and; Rosyada & Hariyanto, 2019)[8, 17]ΛRate of recruitment/birth of human500 × 10^−2^[11, 36]β1Rate of disease transmission from exposed to susceptible[0.45, 0.55][23, 30, 31]β2Rate of disease transmission infected to susceptible[0.45, 0.55][26, 39]μRate of natural death of human[0.01, 0.02][36, 40]δRate of disease induced death0.01[31, 32]αRate of progression of infection from exposed to infected[0.25, 0.35][28, 32]γRate of recovery from infected class[0.07, 0.1][39, 40]γ1Treatment rate to infected class[0.04, 0.06][38, 40]
The birth and death rates in humans (hosts) are relatively low, with a faster population turnover. The transmission rate from the infected class is higher than from the exposed class, i.e., β2ˇ>β1ˇ. Seasonal variations in transmission rates are modeled using sine and cosine functions, capturing seasonal fluctuations and human behavioral patterns. The phase shift between transmission rates influences R0 and significantly impacts the outbreak probability. The value of R0 is computed using the next-generation matrix method and the equations in (25).
In this case, we have selected average parameter values along with a total human population of N(t) = 500. Moreover, the initial numbers of exposed and infected individuals are assumed to be E(0) = 5 and I(0) = 4, respectively. The operation of Brownian motion is illustrated in Fig. 2(a)–(b) using five sample paths. Additionally, these diagrams (a) and (b) provide a detailed view of the exposed and infected hosts. These results offer strong support for the conclusions drawn from the deterministic model.In the first example, the time-periodic transmission rates β1(t) and β2(t) are modeled using sine and cosine periodic functions with the same phase, as shown in equations (26). The other parameter average values are μ = 0.01, α = 0.25, γ1 = 0.05, and δ = 0.01, respectively. Based on these fluctuations, the basic reproduction number R0 is found to lie in the range (World Health Organization https; Ang, 2007). In the non-homogeneous stochastic process, the probability of an outbreak is influenced by the timing of the introduction of the exposed and infected populations, whether at t = 0 or t = 15.
For seasonal effect in all compartments, we also proceed to approach normally distributed parameters β1, β2, α, γ and γ1 [18, 26, 21]. Fig. 5(a) represents a histogram or density plot showing the distribution of contact rates β, β2 for 100 samples. The contact rates follow a normal distribution with a mean of 0.35 and a standard deviation (σ) of 0.04. The plot provides a visual representation of how the contact rates are distributed around the mean value, with the majority of values likely falling within a certain range around 0.35. The confidence intervals are chosen as (μ − 1.5σ, μ + 1.5σ). We observe that within interval [0.3, 0.4], almost 70 samples of β1 and β2 lies.Fig. 5Normally distributed parameters in (a) β1 and β2 with μ = 0.35, σ = 0.04 (b) α with μ = 0.25, σ = 0.04 (c) γ with μ = 0.08, σ = 0.02 and (d) γ1 with μ = 0.05, σ = 0.01.Fig. 5
A histogram or density map of the infection rates for 100 samples is depicted in Fig. 5(b). The infection rates have a mean of 0.25 and a standard deviation (σ) of 0.04, and they are distributed normally. With the majority of values most likely lying within a specific range around 0.25, the plot shows how the infection rates are dispersed around the mean value. We also see that, within interval [0.2, 0.3], almost 70 samples of α lies.
The image is a density plot or histogram that displays the distribution of recovery rates across 100 samples in Fig. 5(c). The recovery rates have a mean of 0.08 and a standard deviation (σ) of 0.02, and they are distributed normally. The plot illustrates the distribution of recovery rates around the mean value, with the majority of values likely falling within a specific range near 0.08. We also observe that over 70 samples of γ lie inside the interval [0.06, 0.1].
By analyzing Fig. 5(d), the illustration shows a histogram or density map of the treatment rates for 100 samples. A normal distribution with a mean of 0.05 and a standard deviation (σ) of 0.01 describes the treatment rates. The plot shows the distribution of treatment rates around the mean value, with the bulk of values most likely falling within a range about 0.05. We also observe that over 70 samples of γ1 lie inside the interval [0.04, 0.06]. These Fig. 5 reflects that, within a certain range, the parameter values oscillates for changing seasonal behaviour.
Because a certain process exhibits recurrent cycles or patterns, seasonal effects with periodic waves in parameters result. These patterns frequently correlate with the cyclical nature of the seasons or other predictable time periods. The periodic waves in the parameters show how a particular quantity or phenomena changes over time, with each wave cycle denoting a distinct interval of time. Numerous natural and artificial processes, including biological rhythms, climate patterns, economic cycles, and social behaviour, exhibit these seasonal influences. To better understand and predict the behaviour of the event under study, we can model these recurring patterns by include periodic waves in the parameters.
amplitude = 0.2, ϕ = 0.25, period = 0.1, mean = 0.25 (d) γ with amplitude = 0.2, ϕ = 0.25, period = 0.1, mean = 0.08 and (e) γ1 with amplitude = 0.15, ϕ = 0.25, period = 0.1, mean = 0.05.
Fig. 6(a) and (b) represents the contact rate with exposed and infected, denoted as β1, β2, respectively, plotted against periodic waves. The parameters of the waves are as the mean value is β1 = 0.35, β2 = 0.4, the amplitude is 0.2, the phase shift is 0.25, and the period is 0.1. In simpler terms, Fig. 6(a) and (b) shows how the contact rates β1, β2 vacillates over time in a cyclical pattern. The waves indicate that the contact rate fluctuates periodically around an average value of 0.35 and 0.4, respectively. The amplitude of the waves is 0.2, meaning the contact rate deviates from the average by up to 0.2 in both positive and negative directions. The phase shift of 0.25 determines the starting point of the wave cycle. Finally, the period of 0.1 indicates the time it takes for one complete cycle of the wave pattern. We observe that, within week 1 to 4, the peak level of β1 reaches 0.42, and the minimum level is 0.28; while the peak level of β2 approaches 0.48, and the minimum level strikes at 0.32.Fig. 6The effect of seasonal succession (periodic transmission rates) in model parameters where, (a) β1 with amplitude = 0.2, ϕ = 0.25, period = 0.1, mean = 0.35 (b) β2 with amplitude = 0.2, ϕ = 0.25, period = 0.1, mean = 0.35 (c) α with.Fig. 6
Fig. 6(c) depicts the infection rate, α, plotted against cyclical waves. The waves have the following a mean value of 0.25, an amplitude of 0.2, a phase shift of 0.25, and a period of 0.1. In plainer terms, the figure depicts the cyclical variation in infection rate across time. The waves show a periodic fluctuation of the infection rate around an average value of 0.25. Since the waves have an amplitude of 0.2, the infection rate can differ from the average in both positive and negative directions by up to 0.2. The beginning of the wave cycle is determined by the 0.25 phase shift. The duration of the wave pattern's whole cycle is shown by the period, which is 0.1.
Fig. 6(d) illustrates the connection between periodic waves and the recovery rate γ. The waves have the following a mean value of 0.08, an amplitude of 0.2, a phase shift (or starting point) of 0.25, and a period (or time between cycles) of 0.1. The graphic shows how the recovery rate varies over time in accordance with the wave pattern by charting the recovery rate gamma with these periodic waves. The precise properties of the waves, such as their peaks, troughs, and oscillations, will dictate the variations in the recovery rate that correspond to them.
Consequently, Fig. 6(e) describes how treatment rate γ1 and periodic waves interact. The wave's average value is 0.05, its amplitude is 0.15, its phase shift (or starting point) is 0.25, and its period (or length of each cycle) is 0.1. The graphic displays how the treatment rate varies over time according to the wave pattern by charting the treatment rate gamma with these periodic waves. The peaks, troughs, and oscillations of the wave's precise shape and characteristics, as well as the variations in the wave's treatment rate, will be determined by these factors.
A Monte Carlo method is utilized to simulate the entire non-linear process, with a time step Δt chosen to be sufficiently small. This ensures that only one of the 13 events, as listed in Table 2, occurs during each time step.
For ∑(t)Δt < 1, event occurrences follow the probabilities in Table 2. For example, the infection probability from the exposed compartment is β1ESΔt, and the probability of no change is 1 − ∑(t)Δt. To estimate the outbreak probability in the non-homogeneous stochastic process, 1000 sample paths are simulated with initial conditions (E(t), I(t)) = (E(0), I(0)), over t ∈ [0, p]. Each path progresses until T > t, where either E(T) + I(T) = 0 or E(T) + I(T) = 50. If the total infected population reaches 50, an outbreak is considered. The proportion q of sample paths reaching zero estimates the extinction probability, while 1 − q estimates the outbreak probability. Larger or smaller threshold values (e.g., 50) are also considered, especially for larger populations and when R0 is not close to 1.
The outbreak probability is compared to a time-homogeneous Markov process, where β1(t) and β2(t) are replaced by their time-averaged values. In this case, Formula (24) is used. Table 4 summarizes the average basic reproduction number R0~ and outbreak probabilities Pˇoutbreak for the time-homogeneous Markov process, considering human death rate μ, disease-related death rate δ, recovery rate γ, and treatment rate γ1.Table 4The average R0 and for the time-homogeneous Markov process the probabilities of an outbreak are computed, that depends on the initial number of exposed and infected humans, (e, i) = (1, 0), (0, 1), or (2, 2) but not on the introduction of time. The average transmission parameters are β1~=0.35, β2~=0.35, μ = 0.01, α = 0.25, γ1 = 0.05, and δ = 0.01.Table 4ParameterPoutbreakβ1β2γR0~**1 − qE1 − qI1 − (qE)^2^(qI)^2^**0.250.250.052.00.74150.76740.99610.26670.26670.05222.110.76660.80040.99730.28330.28330.05442.220.83070.79780.99830.30.310.05672.330.83130.84600.99860.31670.31670.05892.440.82360.83480.99910.33330.34660.06112.55560.85660.82490.99930.350.350.06332.6670.85150.82080.99950.36670.36670.06562.7780.88100.86400.99970.38330.38330.06782.8890.90310.89190.99960.40.40.073.000.85940.89380.9998
Fig. 7(a)–(b) show the impact of seasonality on the model dynamics, incorporating periodic transmission rates β1(t) and β2(t) as in (26). These rates follow sinusoidal and cosine functions, with a period of p=1007 (weekly) and amplitude σ1 = σ2 = 2.5. Subplots in Fig. 7(a)–(b) illustrate the periodic dynamics of the susceptible, exposed, infected, recovered, and treated populations.Fig. 7The effect of seasonal succession in model compartments where, (a) Seasonality's impact in susceptible, infected and recovered population (b) Seasonality's impact in exposed and treated population, where amplitude = 2.5, period = 100/7, ϕ (phase) = 15, β1 = 0.35, β2 = 0.35, γ = 0.08,γ1 = 0.05, α = 0.15 and μ = 0.01.Fig. 7
Fig. 7(a)–(b) display the dynamics of the SEIRT human subpopulations, specifically focusing on the exposed and infected populations under seasonality, where the transmission rates from exposed to susceptible and from infected to susceptible are both modeled as sine periodic functions. Initially, the maximum number of infected individuals reaches its peak at around 40 weeks, with approximately 160-170 infected hosts. Simultaneously, the minimum number of susceptible individuals is observed at levels between 20 and 30 within the 30 to 40-week range. As shown in Fig. 7(a), the peak of the infected human population wave is notably higher than that of the exposed human population. As in Fig. 7(b), we observe that, the exposed class strikes peak level in time range 0 to 30 weeks. The treatment class gradually increase by oscillating up-to 0-50 weeks, then fluctuate occurs around 200-210 level. Consequently, same scenario occurs in recovered class also, where gradual increase occurs within 0 to 50 weeks. After that, when exposed and infected subpopulation moves under control, the recovered class fluctuates around 300-320 level.
Additionally, similar trends can be observed in Fig. 7(a)–(b), where the transmission rates β1(t) and β2(t) are modeled as time-dependent cosine and sine functions, respectively. In these graphs, the numbers of susceptible, exposed, and infected individuals remain relatively stable between 30 and 50 weeks, as shown in Fig. 7(a)–(b). The peak values of the infected and exposed populations are reached after approximately 30 weeks, influenced by the periodic transmission rates in (26) as well as the effects on other parameters, such as α, γ, and γ1. Moreover, the peak number of infected and exposed human population are calculated as around 150-180 in time interval 30-50 weeks approximately, as a consequence of (26).
Furthermore, in Fig. 8(a)–(b), we observe the impact of the modifications in (26) on the parameters β1, β2, α, γ, and γ1. In this scenario, the transmission rates from the exposed and infected classes are modeled as sine and cosine seasonal functions, while the recovery rate γ and treatment rate γ1 follow cosine periodic functions with a period of p=1005. As a result of the periodic changes in (26), the maximum number of infected individuals reaches approximately 280 in week 50, while the peak number of exposed individuals is observed to be around 250 in week 33. According to the average parameter values from Table 3, this is the only case where the numbers of exposed, infected, and susceptible subpopulations fluctuate within higher ranges throughout the simulation. Lastly, Fig. 8(a)–(b) illustrates the seasonal dynamics of the model (1), incorporating infection rate α and recovery rate γ, both modeled as sine functions. In this instance, the amplitude for each seasonal parameter is set to 2.2, with a phase shift of ϕ = 10. This scenario predicts a higher number of exposed and infected individuals between weeks 30 and 40, which occurs relatively faster than in other cases. By controlling with recovery and treatment rates, the exposed and infected subpopulation gradually decrease and varies in lower range. Also, the susceptible, recovered, treatment population vacillates in 200-250, which satisfactory to control the outbreak.Fig. 8The effect of seasonal succession in model compartments where, (a) Seasonality's impact in susceptible, infected and recovered population (b) Seasonality's impact in exposed and treated population, where amplitude = 2.2, period = 100/5, ϕ (phase) = 10, β1 = 0.35, β2 = 0.35, γ = 0.08,γ1 = 0.05, α = 0.16 and μ = 0.02.Fig. 8
Fig. 9(a)–(c) depicts different scenarios of probability of an outbreak in exposed and infected compartments by oscillation of R0 along with β1, β2, γ1. As for Poutbreak(0)=0.7415,Poutbreak(5)=0.8236.The influence of transmission rates β1(t) and β2(t) on the probability of an outbreak is demonstrated in Fig. 9(a)–(c). In Fig. 9(a) and (b), the transmission rates β1(t) and β2(t) are nearly in phase, as described by equation (26). In contrast, Fig. 6 highlights a phase shift between the transmission rates β1(t) and β2(t), along with time-varying parameters α(t), γ(t), and γ1(t). In Fig. 9, the periodic probability of an outbreak Poutbreak for the exposed population is shown in the left column, while the periodic probability of an outbreak Poutbreak for the infected population is shown in the right column. The branching process estimates are computed and plotted for each time step, dividing the time into intervals t = 0, 1, 2, …, 10. In Fig. 9(a) and (b), the blue curves represent the probability of an outbreak for the initial exposed and initial infected populations in scenario 1, while the magenta curves depict the probability of an outbreak in scenario 2. Red markers are used to indicate variations. It is observed that the scenario 1 curves are smooth, whereas scenario 2 shows fluctuations in the probability of an outbreak, reflecting the variations in the density of the exposed and infected compartments.Fig. 9The outbreak probabilities Poutbreak(t) are plotted for initial conditions of either (E, I) = (1, 0) or (0, 1), with the basic reproduction number R0∈[2,3.5]. (a) Probability of transmission from the exposed to susceptible compartment, (b) Probability of transmission from the infected to susceptible compartment, and (c) Total transmission probability from the infection compartments (E, I) to susceptible, considering two scenarios.Fig. 9
In Fig. 9, the shape changes of Poutbreak(t) are clearly visible in relation to the two transmission rates. The extrema of the outbreak probability are aligned with the extrema of the transmission rates. The maximum probability of an outbreak, which is 0.9031 and 0.8919 for the exposed and infected classes respectively, occurs when β1 = β2 = 0.3833, with R0=2.889. The shift in the peak points (either to the right or left) is mainly influenced by whether the outbreak is initiated by the exposed or infected host. When the outbreak is initiated by an infected host, the shift to the left is determined by the transmission rate β2, while if an exposed host initiates the outbreak, the leftward shift depends on the transmission rate β1. Fig. 9(c) expresses, the scenario 1 and 2, where the total average of probability of an outbreak (combination qE and qI) is displayed. More importantly, in scenario 2 more fluctuation occurs for variational density in exposed and infected subclass.
The approximation indicates that the likelihood of a disease outbreak follows a periodic trend, which is affected by the timing of introducing either an exposed or infected host into the susceptible population. The three numerical examples, each with two scenarios and periodic transmission rates β1, β2, demonstrate that the probability of an outbreak is highest at certain times. These results also highlight the relationship between these times and the peak transmission rates for both exposed and infected individuals. From Fig. 9, we observe that when the infected host is introduced, the probability of an influenza outbreak is significantly higher. In contrast, when the transmission starts from an exposed individual, the outbreak probability remains lower compared to the case where the infected host initiates the disease spread.
The Lévy jump model introduces stochastic dynamics with sudden, unpredictable changes, capturing real-world phenomena like epidemics with abrupt shifts. In the model, this approach integrates jumps into the standard compartmental framework, enriching its ability to represent complex disease dynamics (Kiouach and Sabbar, 2019, 2022; Zhou & Zhang, 2016). By formulating the Lévy jump-driven stochastic differential equation (SDE), we aim to analyze and simulate the model both analytically and numerically. The dynamics of the SEIRT model with stochastic effects can be expressed (6.1)dS=Λ−(β1E+β2I)S−μSdt+σSSdWS+JSdNS,dE=(β1E+β2I)S−(α+μ)Edt+σEEdWE+JEdNE,dI=αE−(μ+δ+γ+γ1)Idt+σIIdWI+JIdNI,dR=γI−μRdt+σRRdWR+JRdNR,dT=γ1I−μTdt+σTTdWT+JTdNT.Here, dWX represents the Wiener process component for compartment X, introducing continuous stochastic fluctuations. dNX is a Lévy jump term representing discontinuous changes (jumps) in the system with jump intensity (rate of occurrence) λX. Further, JX is the jump size for compartment X, representing abrupt changes due to external shocks (e.g., policy changes, sudden outbreaks), which can follow a specified distribution (e.g., normal, exponential). σX denotes magnitude of the Gaussian noise affecting compartment X. Moreover, NX represents Poisson counting process for jumps in compartment X.
Theorem 6 (Global Stability of the Equilibrium Point) Consider the Lévy jump model described by the stochastic differential equations explained in (27).
Suppose the system admits a unique equilibrium point E∗ = (S∗, E∗, I∗, R∗, T∗) that satisfies the deterministic equilibrium conditions.
If the following conditions 1.The basic reproduction number R0 < 1, where R0 is defined as the spectral radius of the next-generation matrix for the deterministic system.2.There exists a Lyapunov function V(X) for X = (S, E, I, R, T) LV(X)=∂V∂t+∇V⋅f(X)+12TrG⊤∇2VG+∫RnV(X+J)−V(X)−∇V⋅Jν(dJ)≤0,where f(X) represents the drift term, G is the diffusion matrix, ν(dJ) is the Lévy measure, and J denotes the jump size.3.The Lyapunov function V(X) is radially unbounded and positive definite, i.e., V(X) → ∞ as ‖X‖ → ∞.
Therefore, the equilibrium point E∗ is globally asymptotically stable in terms of probability.ProofWe will prove the global asymptotic stability of the equilibrium point E∗ = (S∗, E∗, I∗, R∗, T∗) for the stochastic SEIRT model step-by-step. The stochastic SEIRT model with Lévy jumps is given by (27). The equilibrium point E∗ = (S∗, E∗, I∗, R∗, T∗) f(E∗)=0,where f(X) represents the drift term of the model. Now, let us define a Lyapunov V(X)=12(S−S∗)2+(E−E∗)2+(I−I∗)2+(R−R∗)2+(T−T∗)2,where X = (S, E, I, R, T). Clearly, from the •V(X) ≥ 0 for all X, and V(X) = 0 if and only if X = E∗.•V(X) is continuously differentiable.
The stochastic generator L acting on V(X) is given LV(X)=∇V⋅f(X)+12TrG⊤∇2VG+∫RnV(X+J)−V(X)−∇V⋅Jν(dJ),where, ∇V is the gradient of V(X). Meanwhile, f(X) is the deterministic drift vector. Further, G is the diffusion matrix. ν(dJ) is the Lévy measure for jumps. Moreover, J is the jump size. We will compute each term separately.
The drift term ∇V⋅f(X)=∑i=15∂V∂Xifi(X).From the definition of V(X), we compute,∂V∂Xi=Xi−Xi∗,for Xi∈{S,E,I,R,T}.Thus we obtain,∇V⋅f(X)=∑i=15(Xi−Xi∗)fi(X),where fi(X) is the drift term for each compartment. For example, for S, the expression becomes,fS(X)=Λ−(β1E+β2I)S−μS.Similarly, we substitute,(28)fE(X)=(β1E+β2I)S−(α+μ)EfI(X)=αE−(μ+δ+γ+γ1)IfR(X)=γI−μRfT(X)=γ1I−μTAfter expanding each term we have,(S−S∗)fS(X)=(S−S∗)Λ−(β1E+β2I)S−μS.Similarly, we can compute,(29)(E−E∗)fS(X)=(E−E∗)β1E+β2IS−(α+μ)E.(I−I∗)fS(X)=(I−I∗)αE−(μ+δ+γ+γ1)I.(R−R∗)fS(X)=(I−I∗)γI−μR.(T−T∗)fS(X)=(I−I∗)γ1I−μT.Simplify similarly for fE, fI, fR, and fT. Therefore, substitute back into ∇V ⋅f(X).
Now, the diffusion term is,12TrG⊤∇2VG.Since V(X) is quadratic, we have,∇2V=diag(1,1,1,1,1).Now, the trace simplifies to,12TrG⊤∇2VG=12‖G‖2,where ‖G‖2=∑i=15σi2Xi2, with σi being the diffusion coefficients.
The jump contribution term is,∫RnV(X+J)−V(X)−∇V⋅Jν(dJ).Using a second-order Taylor expansion for small jumps,V(X+J)≈V(X)+∇V⋅J+12J⊤∇2VJ.By substituting, we obtain that,V(X+J)−V(X)−∇V⋅J=12J⊤∇2VJ.Thus, we have,∫RnV(X+J)−V(X)−∇V⋅Jν(dJ)=∫Rn12J⊤∇2VJν(dJ).
Next, combining all the terms,LV(X)=∇V⋅f(X)+12∑i=15σi2Xi2+∫Rn12J⊤∇2VJν(dJ).For global stability, we have,LV(X)≤0∀X.Since f(X), σi2, and ν(dJ) are bounded, X converges to E∗ (Imkeller & Pavlyukevich, 2006; Teuerle & Jurlewicz, 2009). Meanwhile, LV(X)≤0 for all X, V(X) decreases along trajectories, ensuring that X converges to E∗. □
Fig. 10 (a) depicts, S(t) class decreasing due to infections caused by exposed (E(t)) and infected (I(t)) individuals at rates β1 and β2, respectively. The term μS and Λ accounts for gradual decrease and increase scenario with randomness. Meanwhile, the I(t) compartment gradually increasing through progression from E(t) at rate α. Decreases occur via recovery (γ), progression to severe illness (γ1), or mortality (δ + μ). Further, the R(t) compartment increasing at a rate γ due to recovery from infected individuals (I(t)). The population in R(t) decreases due to natural mortality at rate μ. The drift term governs gradual deterministic changes, while σSSdWS, σIIdWI and σRRdWR introduces random diffusion (Stochastic noise term) in infection progression. Jump intensity JS, JI, JR and jump factor dNS, dNI, dNR model sudden external perturbations, such as population shifts or large outbreaks.
Fig. 10(b) illustrates the E(t) compartment increasing due to interactions between S(t) and infectious states (I(t), E(t)) at rates β1 and β2. The exposed population decreases due to progression to infection at rate α and natural mortality at rate μ. Moreover, the T(t) compartment tracks individuals who have progressed to a more severe condition or are transferred to another category, increasing at a rate γ1 due to progression from infected individuals (I(t)). The population in T(t) gradually decreases for the effect of mortality at rate μ.Fig. 10The curves display five sample paths illustrating the impact of Lévy jumps in the SEIRT model, where (a) shows the dynamics of susceptible, infected, and recovered populations, and (b) shows the dynamics of exposed and treated populations. The average parameter values are sourced from Table 3, with noise σ = 0.04, jump intensity λX = 0.025, and jump size JX = 0.02.Fig. 10
Stochastic noise is captured by σEEdWE, σTTdWT, while jumps (JEdNE, JTdNT) model sudden escalations, such as emergency transfers or unforeseen transitions. Drift accounts for the typical progression, and diffusion reflects the inherent randomness, with jumps reflecting abrupt changes. Both these Fig. 10 (a) and (b) reflects, E(t) and I(t) population increases to pick level after 23 weeks, after that decreases. After 35 weeks all compartments converges to a stable level.
The comparison scenario of Lévy jumps with Brownian Motion is reflected in Fig. 11. In Fig. 11(a), (b) and (c), the S(t), E(t), and I(t) compartment shows distinct patterns under Brownian motion, Lévy jumps, and deterministic models. In the Brownian Motion case, the addition of noise σS, σE, σI causes small, continuous fluctuations around the deterministic trajectory, reflecting environmental or demographic variability. On the other hand, Lévy jumps, characterized by jump intensity λS, λE, λI, jump size JS, JE, JI, and jump factor dNS, dNE, dNI, introduce abrupt declines or increases, modeling external outbreaks. These jumps make the Lévy model more biologically realistic. While Brownian Motion maintains smoother trajectories with convergence to deterministic curves, Lévy jumps highlight real-world unpredictability, deviating from the deterministic model but capturing essential transient dynamics. The deterministic curve lacks the variability of Brownian Motion or the more realism of Lévy jumps, making the latter more suitable for capturing episodic and rapid changes in the S(t), E(t) and I(t) population.Fig. 11The comparison of Lévy jumps effect with Brownian Motion in the SEIRT deterministic curve where (a) behaviour of susceptible, (b) behaviour of exposed, (c) behaviour of infected (d) behaviour of recovered and (e) behaviour of treated population. The average parameter values are taken from Table 3, with noise σ = 0.04, jump intensity λX = 0.025, and jump size JX = 0.021.Fig. 11
Fig. 11(d), (e) illustrates that, under Brownian motion, the R(t) compartment exhibits smooth fluctuations due to recovery (γ) and natural mortality (μ), with noise σR adding randomness. Lévy jumps, with jump intensity (λR) and jump size (JR), model sudden spikes or declines in recovery, reflecting abrupt external factors like vaccination campaigns. The deterministic model provides a baseline for steady recovery trends but lacks the realism captured by Brownian Motion or Lévy jumps. Similarly, in Brownian Motion case, the T(t) compartment shows gradual changes influenced by transitions (γ1) and mortality (μ), with σT introducing stochastic variability. While in Lévy jumps, defined by (λT) and (JT), capture sudden increases in severe cases due to external shocks, such as rapid disease progression. While deterministic trajectories highlight average trends, Brownian Morion and Lévy jumps implements better capture the dynamic and unpredictable nature of severe disease transfers.
The stochastic version of the SEIRT model can be described by the following SDE:dS(t)=Λdt−(β1E(t)+β2I(t))S(t)dt−μS(t)dt+σdW1(t)dE(t)=(β1E(t)+β2I(t))S(t)dt−(α+μ)E(t)dt+σdW2(t)dI(t)=αE(t)dt−(μ+δ+γ+γ1)I(t)dt+σdW3(t)dR(t)=γI(t)dt−μR(t)dt+σdW4(t)dT(t)=γ1I(t)dt−μR(t)dt+σdW5(t)where, the state variables and parameters are described in Table 1
σ is the noise intensity and dW1, dW2, dW3, dW4 and dW5 are independent Weiner processes representing the stochastic noises.
Our goal is to fit the SDE model described in Section 3 to the real data. We have observed data points for I(t) over time. We want to find the parameter values β1, β2, γ, γ1, μ, δ and noise term σ that best fit the data. To achieve this, we apply the least squares method to minimize the sum of the squared differences between the model predictions and the observed minβ1,β2,α,γ,γ1,μ,δ=∑iIiobs−I(ti;β1,β2,α,γ,γ1,μ,δ)2Here, Iiobs are the observed data points at time at time ti and I(ti; β1, β2, α, γ, γ1, μ, δ) are the corresponding model predictions. Estimating the noise term σ involves finding the value that best captures the variability between the model and the real-world data. This can be done alongside the parameter estimation using nonlinear least squares optimization. The process of estimating noise σ involves iteratively adjusting its value while optimizing other parameters to minimize the difference between model predictions and observed data. Estimating σ accurately can be challenging, as it reflects the level of randomness in the system and is influenced by factors such as measurement errors and unaccounted-for sources of variability. To estimate the noise intensity, we calculated the sample variance of the differences between the simulated and observed trajectories. We also incorporated the estimated noise term into the SDEs and refine the parameter estimates based on the noise-adjusted model.
To examine the current trend using our proposed model, we analyzed recent influenza infection data from Mexico, covering the period from October 1, 2020, to March 31, 2023. The data was sourced from the CDC and WHO websites (Centers for Disease Control and; World Health Organization https). We considered a total of 120 weekly data points to evaluate the model's applicability and project future disease trends. The initial population values were set as S(0) = 990, E(0) = 5, I(0) = 4, R(0) = 1, and T(0) = 0. The predicted outcomes of the model, based on this real data, are shown in Fig. 12(a) and (b), illustrating the weekly reported cases. Here, the noise term σ is estimated for the exposed and infected compartments, while σ1 is estimated for the susceptible, recovered, and treatment classes using the nonlinear least squares method. A linear regression analysis of the average parameter values in Table 3 shows that the model aligns remarkably well with the weekly reported cases of infection, indicating a strong fit. Consequently, the model demonstrates an accuracy of approximately 76% in tracking the original confirmed cases. Mexico, located in North America and bordering the United States, is a key region of interest. According to a CDC report, during the 2018-2019 influenza season, Mexico saw nearly 7500 reported cases. By March 13, 2021, however, only three cases had been recorded for the 2020-2021 season. In the depicted period, the number of influenza cases reached its peak during the 2015-2016 season, with nearly ten thousand reported occurrences. Mexico, the country most severely affected, has documented 176 fatalities due to the new strain of the H1N1 virus. Cases have been observed globally, primarily among travelers from Mexico, though they have generally been mild, and individuals have recovered through appropriate vaccination and treatment. By the beginning of 2021, reported cases of influenza in Mexico stood at 5.2%, with a death rate of 1.2%. In response, authorities initiated a vaccination program that is still ongoing in hospitals, community vaccination sites, clinics, and pharmacies. After receiving the recommended doses, the immunizations proved to be effective.Fig. 12Model data fitting to Mexico with weekly infected cases, where the fitted parameter sets β1 = 0.35, β2 = 0.35, μ = 0.2, δ = 0.01, γ = 0.05, α = 0.25 and Noise σ = [0.045, 0.5], σ1 = [0.01, 0.03].Fig. 12
Afluria Quadrivalent, Fluarix Quadrivalent, FluLaval Quadrivalent, and Fluzone Quadrivalent are additional approved vaccines. These vaccines are safe for administration to infants as young as six months. The efficacy of these vaccines is approximately 66%, 64%, 67%, and 62%, respectively. In Fig. 12(a) and (b), two waves of weekly reported infected data from different regions of Mexico are depicted. The figures illustrate the correlation between the model's predictions and the actual reported data. The results indicate that the proposed model is applicable. When analyzing the total number of cases per 1000 population, Fig. 12(a) and (b) show that seasonal flu cases in Mexico are recorded sporadically each week, ranging from 0 to 20 weeks. By week 20, around 500 cases were reported. The infection gradually increased from week 21 to week 40, with the highest number of reported cases (720 per 1000) occurring during weeks 33 and 37. From weeks 30 to 50, the disease conditions were particularly severe. Thanks to proper care and isolation measures, the epidemic began to subside after week 50. After 80 weeks, the reported cases per thousand ranged from 100 to 150. Notably, when the noise terms are set to σ = 0.05 and σ1 = 0.02, the model's data aligns satisfactorily with the real-world data.
Fig. 12(a) and (b) illustrate that insightful conclusions were drawn from our model's analysis of influenza cases in Mexico. The data provides a comprehensive overview of the influenza situation and offers valuable insights into the virus's spread and impact. By fitting the data to our model, we were able to make accurate predictions and identify key patterns and trends. These results can assist researchers, policymakers, and public health experts in formulating effective strategies to combat influenza and mitigate its effects on the population. The findings of our model present a compelling and informative account of influenza data in Mexico, contributing to ongoing public health initiatives and disease surveillance efforts.
In this study, we applied the SEIRT epidemic model to represent host populations in the influenza model. By taking into account all key aspects of the disease progression, we have developed both stochastic and deterministic dynamics for the influenza model. Our Brownian motion, Lévy jump stochastic model allowed us to develop a number of theoretical conclusions and graphically depict them in terms of the SDE. Effectively, there is a satisfactory level of agreement discovered between the deterministic and stochastic models' dynamic behavior. It is also demonstrated that, after a certain period, the number of exposed and infected individuals will gradually decrease. Thus, the effectiveness of the strategy to reduce and ultimately eliminate disease outbreaks is ensured by our model, which incorporates both deterministic and stochastic elements. In this study, we integrated Brownian motion, Lévy jump and seasonality into an Influenza SEIRT model to enhance the representation of disease dynamics. Our simulations revealed that these additions significantly impact the model's ability to mimic real-world influenza dynamics. Notably, the model successfully captured sporadic outbreaks, seasonal peaks, and the efficacy of treatment interventions. By incorporating a noise term and jump intensity, we effectively captured the inherent stochasticity in disease transmission and calibrated the model parameters to observed data. The results demonstrated a notable improvement in the model's ability to reproduce real-world influenza dynamics. The data fitting approach using Mexico data allowed us to refine the model, aligning it closely with empirical observations of infection spread and treatment outcomes. The comparison results of each scenario highlight the significance of incorporating both stochastic and seasonal factors in influenza modeling, leading to a deeper understanding of disease transmission and offering important insights for crafting effective public health interventions and response strategies. Additionally, after analyzing various theorems on the extinction and persistence of seasonal flu, it is confirmed that the treatment rate can be a critical factor in reducing the peak infection level.
For a non-homogeneous, time-periodic stochastic model, we showed that the backward Kolmogorov differential equations for the non-homogeneous multi-type branching process approximation can be used to derive a system of differential equations for estimating the asymptotic probability of an outbreak, Poutbreak(t) for t ∈ [0, p] (see (23), (24)). The number and timing of introductions of infected or exposed hosts into a susceptible population, along with periodic transmission rates, determine the periodic outbreak probability. Analytical and simulation results reveal that the peak transmission period does not align with the peak outbreak risk; instead, outbreak probability peaks earlier than transmission rates.
The approach assumes large population sizes and a small number of initial infections, with independence among initial infections. When R0<1, extinction is guaranteed (subcritical case). For R0>1, the process is supercritical, leading either to extinction or unbounded growth. A sufficiently large host population (e.g., N(t) = 1000) is required for the branching process approximation to be valid; populations below 100 may yield less accurate estimates.
These results are broadly applicable to non-homogeneous stochastic epidemic models with periodic transmission, relevant for diseases like malaria, dengue, Zika, Lyme disease, chikungunya, leishmaniasis, and influenza, which are influenced by seasonal factors. Estimating periodic transmission rates by understanding host-pathogen responses to environmental signals is crucial for predicting outbreak risks. Public health interventions must account for seasonal effects and climate change, and the models discussed here offer valuable tools for designing control strategies, preventive measures, and treatment policies.
Kazi Mehedi Mohammad: Writing – original draft, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. Taufiquar Khan: Writing – review & editing, Validation, Software, Resources, Investigation, Funding acquisition. Md Kamrujjaman: Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Formal analysis, Conceptualization.
Not applicable.
The datasets generated and/or analyzed during the current study are available in the Md Kamrujjaman [jamanmd] repository, https://github.com/jamanmd/Influenza.
No consent is required to publish this manuscript.
Not applicable.
The project was partially supported by the University Grants Commission of Bangladesh.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.