A stochastic SIQR epidemic model with Lévy jumps and three-time delays.

Isolation and vaccination are the two most effective measures in protecting the public from the spread of illness. The SIQR model with vaccination is widely used to investigate the dynamics of an infectious disease at population level having the compartments: susceptible, infectious, quarantined and recovered. The paper mainly aims to extend the deterministic model to a stochastic SQIR case with Lévy jumps and three-time delays, which is more suitable for modeling complex and instable environment. The existence and uniqueness of the global positive solution are obtained by using the Lyapunov method. The dynamic properties of stochastic solution are studied around the disease-free and endemic equilibria of the deterministic model. Our results reveal that stochastic perturbation affect the asymptotic properties of the model. Numerical simulation shows the effects of interested parameters of theoretical results, including quarantine, vaccination and jump parameters. Finally, we apply both the stochastic and deterministic models to analyze the outbreak of mutant COVID-19 epidemic in Gansu Province, China.

Introduction

Infectious disease is a strong antagonist, posing a threat to human life and health. Recently, the outbreak of the novel coronavirus (COVID-19) had a huge impact on social life all over the world. The control and eradication of the disease have been an urgent problem and received considerable attention from the researchers. Mathematical modeling is one of the most commonly used methods for analyzing the transmission mechanisms as well as controlling of various epidemic diseases. Kermack and Mckendrick [1] used a rigorous mathematical approach and for the first time, they formulated a susceptible-infected-recovered (SIR) epidemic model to discuss the spread of contagious diseases. Following that, numerous researchers used ordinary differential equations for modelling different types of infectious diseases, for instance, one can see [2], [3], [4]. During the transmission of many diseases, quarantine and vaccination are the two most effective methods to prevent the epidemic from spreading further, such as COVID-19. Hence, more and more compartments and factors that inhibit the spread of diseases were considered in the SIR model. Suppose that the total host population is partitioned into susceptible, infectious, quarantined and recovered individuals, respectively, whose size is , , , and at any time . Ma et al. [5] formulated and analyzed an SIQR model assuming the hybrid strategies and the bilinear incidence rate in the following formwhere the parameter is the constant recruitment into the susceptible compartment, is the rate at which the disease is spreading in the population, is the vaccination rate of the susceptible individuals, represents the elimination rate of the infective class and the natural death rate is denoted by which is constant for all compartments. The parameter denotes the recovery rate of individuals that remain infectious for a period of , the notations and are the respective disease-related death rates in the infected and quarantined populations. The rate at which the infected population becoming quarantine is denoted by , whereas, is the rate at which quarantine population is getting recovery. For biological purposes, we assume that all of the parameters used in the model are constants and positive.

In practice, the history related to an epidemic can play a significant role in studying the dynamics of certain infectious diseases. That is because, at any time , the dynamics of diseases are subject to some previous conditions. Time delay is the most generally used tool to simulate the infectious period in the models which describe dynamics of epidemics. Such conditions occur frequently and extensively in the actual world such as temporary immunity and incubation period. Therefore, it is natural to introduce time delay into the SIQR model describing the dynamics of an infectious disease [6], [7], [8]. Lu et al. [9] proposed an SIQR epidemic model with one time-delay to describe the spreading behavior of COVID. However, there are two limitations of this model: (i) It fails to consider the effect of COVID-19 vaccination. (ii) There exist other time delays in the susceptible and quarantined periods, similar to the infectious period. Inspired by the works of Ma et al. [5], [9], we introduce three-time delays into the SIQR model (1) with vaccination, and establish a novel model aswhere is positive and stands for the length of immunity period of the vaccinated individuals, the term represents the probability that a susceptible individual is vaccinated at time and he/she is still alive at time . Similarly, the second and third delay terms and are positive and these reflect the durations of the immunity period of individuals being recovered from infected and quarantined compartment, respectively. In the same way, the quantities and denote the probabilities of infected and quarantined individuals who are respectively immunized at time and , lose their immunity and still alive at time . The notion is the rate of recovered people who forget to re-vaccinate before losing immunity. Since the first three equations in model (2) do not depend on the fourth equation, thus, we have the following simplified form of the modelFrom model (3), we can easily obtain the basic reproduction numberwhich determines whether the disease is tend to extinct or not. The deterministic model (3) can describe the general trends of epidemic and has a good prediction effect for large-scale infectious diseases. The dynamic properties of this model will be discussed later in this work.

However, epidemics are inevitably affected by the environmental noises, which are the important components to be taken into account by mathematical models. Therefore, it is suitable to use the tools of stochastic modelling to describe such characteristics of infectious diseases, and perform better on local infectious diseases or a small number of infections, comparing to the deterministic model (3). In this regard, various stochastic epidemic models were formulated and analyzed to reveal the influence of environmental noises [10], [11], [12], [13], [14], [15], [16]. El Fatini et al. [17] proposed a triple delayed SIQR model with white noises. Further, Pitchaimani [18] studied the triple delayed SICR model with general incidence rate. In these results, the authors only assumed that the stochastic models were driven by the white noises. For some highly infectious diseases, the environmental noises may lead to severe fluctuations in the number of infected people, such as COVID-19 [19], [20], Tuberculoses [21], etc. Thus, it is inappropriate and unrealistic to prefer the deterministic over stochastic approach while modeling such epidemics. Based on these facts, numerous researchers take into account the stochastic epidemic modeling driven by Lévy jumps [22], [23], [24], [25], [26], [27], [28], [29], [30]. Recently, Koufi et al. [31] proposed a stochastic SIQR epidemic model and discussed the effect of Lévy jumps as well as utilized the incidence rate of Beddington-DeAngelis type. The model of Koufi et al. [31] is indeed a great contribution to stochastic epidemic model from various perspectives. However, the model ignored the effect of time delays which arise naturally. Besides that, they did not investigate the asymptotic properties of solutions around the equilibrium points.

Motivated by the aforementioned work, in this paper, we employ three white noises and a stationary Poisson point process as the driven jump process into model (3). Let us consider a complete probability space , , , with a filtration that is right continuous as well as increasing. Also, contains all -null sets. Suppose that is a Poisson counting measure with characteristic measure on measurable subset of , satisfying . Define for , which are used for reflecting the impacts of random jumps. The compensated random measure is denoted by and is defined by . Based on model (3), we propose a stochastic epidemic model with Lévy jumps given by the following systemwhere denote the left limits of , and , respectively. The notions are the mutually independent standard Brownian motions. Similarly, are the perturbation volatilities.

To the best of our knowledge, due to the complexity of constructing appropriate Lyapunov functions and related calculations, there are very few studies on the asymptotic properties around the equilibrium points of delayed models. The present work is indeed a great novelty for the stochastic delay approach and for the asymptotic properties around the equilibrium points in stochastic models. The rest of this paper is organized as follows. We review some definitions and notation which are used in the later parts of the study in Section 2. In Section 3, we firstly prove that there is a unique positive global solution of the stochastic SIQR model (5). Then, we investigate the dynamic properties of stochastic solutions around the equilibrium points of the model (3). In Section 4, some numerical solutions are presented to validate the obtained analytical findings. In Section 5, the spread of mutant COVID-19 is investigated through our proposed models in Gansu Province, China. A brief conclusion is given in Section 6.

Preliminaries

Consider the following three-dimensional stochastic differential equation with Lévy jumpswith the initial value where . The operator acting on the equation  (6) is defined by

Let be the Banach space of all continuous mappings which is equipped with the norm such that . Denote as the family of all nonnegative functions defined on such that they are continuously twice differentiable in and once in . By applying the operator on a Lyapunov function , it follows thatwhereBy utilizing the Itô’s formula, we obtain

Main results

Before studying the dynamical behavior of model (5), it is perhaps the most significant step to show the existence of a global positive solution of the model. Let . Denote and .

For an initial value, there exists a unique positive solutionof model(5)for all time. The solution will remain in the spacewith unit probability.

Obviously, all the coefficients of model (5) are locally Lipschitz continuous. Hence, the existence of a unique local solution is sure for all and an initial value that lies within the space . Here, the notion represents the duration of the explosion. To show that this solution is global as well, it is sufficient to prove that almost surely. To prove this, we consider a sufficiently large positive integer such that belong to the interval , respectively. Further, for every integer greater or equal than , the stopping time is defined aswhere and is the null set. Surely, is a monotonically increasing function of . Set , which implies almost surely. If we prove that , it will reflects that and ultimately assures that lies within the space a.s. for all . If this statement is false, then there exists a constant and from the interval (0,1) such that

As a result, there must exist an integer such thatTo proceed further, a -function from the space to is assumed aswhere is a positive constant to be determined later. By using the Itô formula (8), we havewhereChoosing and keeping in mind the fact thatwhere is a constant, we haveHence, it yields thatBy integrating both sides of Eq. (11) from 0 to and then assuming the expectation, we have

For any , we have . Then, . Consequently, at least one of the classes equals either or for each . ThusFrom the above, we havewhere is the indicator function of . If , thenThus, we have a.s. and the theorem is proved. □

Theorem 1 reveals the existence and uniqueness of the stochastic solution of model (5). In the model (3), the disease-free equilibrium is denoted byIt is handy to prove that if , the equilibrium is both locally and globally asymptotically stable. In the case of stochastic models like the system (5), researchers have putted a question mark on the existence of . Consequently, one can explain the behavior of solution of the models around the disease-free equilibrium. Particularly, it is interesting to study that what kind of changes will appear around the disease-free equilibrium. The following theorem aims to answer such and related questions. For convenience, denote

Assume thatis a solution of model(5)with the initial value. Ifand, thenwhere

Let in model (5). Thus, model (5) is equivalent to the following equationswhere and are positive functions. Define another -function of the formwhere is a positive real number to be calculated later. By utilizing the formula due to Itô (8), we havewhere the term in equation (13) can be written as Choose , that is, . Thus, for , we haveIf we integrate Eq. (13) over the interval and take the expectation, it follows thatBy using the inequality (14) in the above relation, we haveThereforewhere  □

In Theorem 2, we observed that for and , every solution of model (5) will be swinging in the vicinity of disease-free equilibrium of the model (3). Further, the disturbance range is proportional to the values of ’s and ’s for each . In other words, if one assumes very small white and Lévy noise intensities, then every solution of model (5) will cluster very close to . In particular, if , and , then is also the disease-free equilibrium of model (5)

According to the proof process in Theorem 2, we know that is negative definite and is positive definite. It means that the trivial solution of model (5) is stochastically asymptotically stable.

If,and, then the disease-free equilibriumof model(5)is stochastically asymptotically stable in the large.

Next, we investigate the asymptotic behavior of the stochastic solution of model (5) around the endemic equilibrium of the model (3) asHere again, we introduce the following notations for the sake of simplicity

Assume thatis a solution of model(5)with the initial value. Ifand, thenwhere.

Define a non-negative function of the formwhere is a positive real number to be calculated later. By the Itô’s formula (8), we havewhereFrom the model (3), the endemic equilibrium satisfies the following equationsTherefore, by substituting these relations in Eq. (16), we haveChoose so that . Utilizing the inequalities and , we have

Here again, if we integrate both sides of the Eq. (15) over the interval , and then take the expectations, we obtain the following inequalityDenote . Then  □

Physically, Theorem 3 indicates that if and , all the solutions of model (5) fluctuates around the endemic equilibrium of the model (3). Here again, the disturbance range of the fluctuations is proportional to the values of and . In other words, if the intensities of Lévy noise and white noise become smaller, each solution of system (5) will remain close in the vicinity of .

Numerical simulations

To verify the obtained theoretical results, different settings of parameters together with numerical simulations are proposed in this section. Model (5) consists of three parts: deterministic part with delays, white noises determined by Wiener process, and Lévy jumps. For the general stochastic model with delays, we usually use the Milsteins higher-order method and Euler-Maruyama algorithm [32]. Similar to [20], the jump part is the compound Poisson process defined by where is a Poisson process with mean . The jump size is independent and identically distributed random variable with distribution function . Here, follows the standard normal distribution, that is, . The Lévy jump measure of is given by . In the simulation, we assumed .

Assume that , , , , , , , , , , . The time step size . The initial value of all the compartments in vectorized form is given as . By using these values of the parameters, the conditions of Theorem 2 and are satisfied. Fig. 1 shows a sample solution of model (5) which clearly fluctuates around the disease-free equilibrium proposed by the underlying deterministic model (3).

Fig. 1

(a) The trajectories of in models (3) and (5) with parameters values defined in Example 4.1. Panel (b) and (c) shows the trajectories of and , respectively.

Assume that and other parameters are similar to those of Example 4.1. Through calculation, we have and . Hence, the condition of Theorem 3 holds. From Fig. 2 , each solution of model (5) fluctuates around the endemic equilibrium of model (3).

Fig. 2

(a) The trajectories of in models (3) and (5). Panel (b) and (c) shows the trajectories of and with parameters values defined in Example 4.2.

Assume that , , and consider other values of the parameters in Example 4.1. It can be seen from Fig. 3 that when and , then by increasing the intensity of the noise will decreases the numbers of infectious and isolated individuals.

Fig. 3

(a) describes the trajectories of in stochastic model (5) with or without jumps, which the parameters values are defined in Example 4.3. (b) describes the trajectories of in the stochastic model (5) with or without jumps, compared with model (3). Here, model (5) with jumps has Lévy noise and white noise. Model (5) without jumps only has white noise.

The values of reflect the trajectories of model (3). The values of correspond the trajectories of model (5). (a) describes the trajectories of the solutions of models (3) and (5) with the isolation rate and other parameter values of Example 4.2. Panel (b) and (c) are the trajectories of models (3) and (5) under the cases and , respectively.

Assume that the isolation rate and other parameters are defined in Example 4.2. Fig. 1 reveals the effect of the rate by . On the other hand, we analyzed the sensitivity of the isolation rate according toNote that has the negative sensitivity index for any positive parameters. This indicates that will decrease as we increase the value of the quarantine parameter. Thus, the increase of will help the elimination of the infectious individuals from the population.

Assume that the vaccination rate and other parameters are defined in Example 4.2. Fig. 5 reflects the effect of different vaccination rates in models (3) and (5). The sensitivity of is expressed byBy substituting all the parameters into the above equation, we observe that has the negative sensitivity index. Hence, decreases when the vaccination parameter increases. Similar to the isolated rate , the vaccination rate can help to cut the number of the infectious.

Fig. 5

The values of reflect the trajectories of model (3). The values of correspond the trajectories of model (5). (a) describes the trajectories of the solutions of models (3) and (5) with the vaccination rate and other parameter values of Example 4.2. Panel (b) and (c) are the trajectories of models (3) and (5) under the cases and , respectively.

Mutant COVID - 19 epidemic

Mutant COVID-19 strains are now spreading geographically much faster than at any time in the history. Studying the spread of COVID-19 epidemic transmission can help better control and prevent the epidemic. Recently, China has suffered another hardest-hit by the virus since the beginning of the 2022. In order to study the outbreak, we obtained the number of existing confirmed cases in Gansu Province from March 11 to April 2 from the Municipal Health Commission (http://wsjk.gansu.gov.cn/). According to the National Bureau of Statistics, the total population of Gansu Province has reached 25,019,831 from the seventh national population census. Thus, and .

The unknown parameters are estimated by the least-square method [33] with the initial value . Table 1 lists the estimated values of parameters of models (3) and (5). Take other parameters from Example 4.1. Through calculation, and satisfying the conditions of Theorem 2. Therefore, the mutant COVID - 19 epidemic of Gansu Province will become extinct over time. Based on all values of parameters, Fig. 6 shows curves of the fitted mutant COVID-19 cases in models (3) and (5), compared with the real data from March 11 to April 2 in Gansu Province. Obviously, it can be seen from Fig. 6 that the stochastic model is considered to be better in fitting the real data.

Table 1

The estimated values of model parameters by least-square method.

Parameter	Physical meaning	Estimated value
μ	Natural death rate	0.065
p	Vaccination rate of the susceptible individuals	0.0015
ϕ	Disease transmission rate	1.01e-8
ρ	Recovery rate of infected individuals	0.0098
π	Rate of isolated individuals who are infectious	0.013
κ	Recovery rate of quarantined population	9.4e-4
ς_1	Disease-related death rates in the infected	0.067
ς_2	Disease-related death rates in the quarantined	0.004
m	Loss of immunity rate of vaccinated individuals	0.029

Fig. 6

The trajectories of the fitted mutant COVID-19 cases in models (3) and (5), compared with the reported existing infections of Gansu Province, China. Model (3) corresponds to a blue dotted line. Model (5) corresponds to a red solid line, and the red asterisks are the real data. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Conclusion

The spread of infectious disease is influenced by various factors including recruitment, transmission, recovery, vaccination, and isolated rates. Time delays and Lévy jumps are also very important for modelling the spread of epidemic. In this paper, we have proposed a stochastic SIQR epidemic model with triple delays and Lévy jumps. Through the theoretical analysis of the stochastic model, the main findings are listed as follows: (i) Using the suitable Lyapunov method with time delays, we obtain the existence and uniqueness of global positive solution of model (5) (see Theorem 1). (ii) If and , the stochastic solution of model (5) fluctuates around the disease-free equilibrium of model (3) (see Theorem 2 and Remark 1). Especially, if , then the solution of deterministic model (3) and stochastic model (5) have the same asymptotic behavior (see Corollary 1). (iii) If and , then the stochastic solution varies around the endemic equilibrium of model (3) (see Theorem 3 and Remark 2). Because of the complexity of the stochastic model, all the theoretical results are obtained by constructing Lyapunov function with time delays. Further, they reveal the relationship of models (3) and (5).

The numerical simulations are used to illustrate our analysis results. Under different parameter settings, we discussed and analyzed the effects of interested parameters in models (3) and (5) according to the basic reproductive number and other conditions. We summarize some useful results: (i) the disease becomes extinct if the transmission rate is smaller (Figs. 1 and (2). (ii) The increasing of noises may lead to a decrease the number of infectious and isolated individuals (Fig. 3). (iii) The infected individuals will decrease if the isolation and vaccination rates increase (Figs. 1 and (5). All the results are based on the fixed Lévy jumps and three time delays. Finally, we apply the proposed models to investigate the spread of mutant COVID - 19 epidemic in Gansu Province. The result reflects that the mutant COVID - 19 epidemic eventually go to extinct.

In our work, there are still some interesting issues which needs further research. Because of the complexity in the proposed models (3) and (5), more of the theoretical properties has not been solved such as the existence of unique ergodic stationary distribution. They will be considered in our further work.