Stochastic model of the transmission dynamics of COVID-19 pandemic

In this paper, we formulate an SVITR deterministic model and extend it to a stochastic model by introducing intensity of stochastic factors and Brownian motion. Our basic qualitative analysis of both models includes the positivity of the solution, invariant region, disease-free equilibrium point, basic reproduction number, local and global stability of disease-free equilibrium point, endemic equilibrium point, and sensitivity. We obtain the stochastic reproduction number and local stability by using twice differentiable Itô’s formula. We prove the global stability of the disease-free equilibrium point by using a Lyapunov function. We determine the sensitivity of the effect of each parameter on basic reproduction number of the model by using a normalized sensitivity index formula. On the other hand, we demonstrate numerical simulation results of deterministic and stochastic models of COVID-19 by using Maple 18 and MATLAB software. Our simulation results indicate that reducing the contact between infected and susceptible individuals and improvement of treatment play a vital role in COVID-19 pandemic control.


Introduction
Mathematical modeling is useful in understanding and analyzing the behavior of infectious disease transmission dynamics in humans and animals. Since the first discovery and identification of coronavirus in 1965, three major outbreaks occurred, which were caused by emerging and highly pathogenic coronavirus. In 2003, in Mainland China, the outbreak of Sever Acute Respiratory Syndrome Coronavirus (SARS-CoV) occurred [5,12]; in 2012 the outbreak of Middle East Respiratory Syndrome Coronavirus (MERS-CoV) in Saudi Arabia occurred [16], and the virus MERS-CoV in South Korea happened [3,10]. Currently, Corona Virus Disease 2019 (COVID- 19) is an infectious disease caused by Sever Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). The infection was noticed in Wuhan, China, where by the first case of infection was identified in December 2019 [2,11,13]. The outbreak was declared a public health emergency of international concern on 30 January 2020 by World Health Organization (WHO).
The World Health Organization renamed COVID-19 as strain Severe Acute Respiratory Syndrome Corona Virus 2 (SARS-CoV-2) on 11 February 2020 [4,20]. As of 17 March 2021, more than 120 million cases of COVID-19 infections, with more than 2 million deaths, had been reported globally [20]. According to Center for Diseases Control (CDC) [1], the major symptoms of COVID-19 are fever, cough or sneeze, and breath shortness, which emerges 14 days after the infection occurred. Public health [18] and WHO [20] claimed that the virus that causes COVID-19 spreads mainly when an infected person is in close contact with another person. Small droplets and aerosols containing the virus can spread from an infected person's nose and mouth when they breathe, cough, sneeze, sing, or speak [21]. The virus may also spread via contaminated surfaces. Based on the spreading nature of COVID-19, WHO declared it as a global pandemic on 11 March 2020 [20].
It is known that there is no known medicine to combat the novel coronavirus pandemic yet [14]. WHO made certain standardized recommendations to prevent the spread of COVID-19. They include frequent washing of hands with soap and alcohol-base sanitizer, using face mask, and avoiding close contact with anyone that has a fever and cough [21]. A COVID-19 vaccine is intended to provide acquired immunity against (SARS-CoV-2). In phase III trials, several COVID-19 vaccines have demonstrated efficacy as high as 95% in preventing symptomatic COVID-19 infections. So far, more than 20 million doses of the AstraZeneca vaccine have been administered in Europe, and more than 27 million doses of the Covishield vaccine (AstraZeneca vaccine by Serum Institute of India) have been administered in India. As of 17 April 2021, 890.31 million doses of COVID-19 vaccine have been administered worldwide by official reports from national health agencies [9].
Mathematical models are important tools to gain a big understanding of the ongoing trends for COVID-19. They are also useful for obtaining a basic reproduction number, determining sensitivities to change in parameter values, estimating key parameters from the data that contribute to identifying trends, making general forecasts, and estimating uncertainties [6]. Epidemiological models play a fundamental role in the study of the dynamics of COVID-19. With regard to the studies carried out so far, a few mathematical modeling studies have been done about transmition of the pandemic. For instance, Wu et al. [22] introduced a susceptible-exposed-infectious-recovered (SEIR) model to describe the transmission dynamics and forecasted the national and global spread of the disease. The model is considered simple and does not incorporate the relapse form recovered class to susceptible class, which is very unrealistic. The author did not introduce any treatment and vaccination classes. A similar study was done by Read et al. [15], but with no vaccination and treatment classes. As COVID-19 is a chronic disease, treatment should be considered to make the work more realistic, which is missing in [15]. Tang et al. [17] proposed a deterministic compartmental model incorporating the clinical progression of the disease, the individual epidemiological status, and the intervention measures. They found that the control reproductive number could be as high as 6.47, and those intervention strategies such as intensive contact tracing followed by quarantine and isolation can effectively reduce the control reproduction number and the transmission risk. Imai et al. [8] conducted computational modeling of potential epidemic trajectories to estimate the size of the disease outbreak in Wuhan, with a focus on the human-to-human transmission. Their findings showed that control measures need to block well over 60% of transmission to be effective in containing the outbreak.
Recently, Iboi et al. [7] developed a mathematical model to determine whether or not a hypothetical imperfect vaccine can lead to the elimination of COVID-19 in the United States. Their study indicated that such elimination is feasible, using the hypothetical vaccine with assumed efficacy of 80%, if the vaccine coverage is adequate enough to achieve herd immunity. In particular, the vaccine coverage needed to achieve herd immunity in US is 90%, whereas the computed herd thresholds for the states of New York and the state of Florida are 84% and 85%, respectively. Therefore from the above discussion we can conclude that human-to-human contact is the potential cause of outbreaks of COVID-19. Thus isolation of an infected human can reduce the risk of future spread of COVID-19. Thus, to overcome those limitations, we conducted the current study to develop a stochastic SVITR mathematical model for the transmission dynamics of COVID-19 pandemic by introducing treated and vaccinated classes.
The paper is organized as follows. In Sect. 2, we formulate and describe the mathematical model. In Sect. 3, we qualitatively analyze the model by examining the equilibrium points. Numerical simulations of the model by estimating the parameters are given in Sect. 4, where the sensitivity of the basic reproduction number on the model parameters is also discussed. Last but not least, conclusions and recommendations of the study are given in Sect. 5.

Model description
The model we have already formulated consists of five compartments: susceptible S(t), vaccinated V (t), infected I(t), treated T(t), and recovered R(t) individuals for all times t > 0. The susceptible are those individuals that are not infected by COVID-19 but can infect in future. Infected individuals include individuals that are couched by COVID-19 and are able to transmit the disease to a susceptible one. Treated individuals are those who can be treated from the disease and cannot transmit to other individuals. After treatment class, some individuals are recovered from the disease. The susceptible populations are increased by recruitment rate ϕ, either by birth or immigration. Those individuals can be vaccinated at rate θ and become infected by COVID-19 with contact rate α. The infected individuals are decreased by disease causing death rate τ . The infected individuals can get treatment at rate δ, and those individuals that are in treatment can recover from the disease at rate ρ. The vaccinated and recovered individuals can lose their temporary immunity at rates σ 2 and σ 1 , respectively, becoming susceptible again. We have natural causing death rate μ for the whole population.
Our model is governed by the following assumptions: All parameters are nonnegative, the total population size is constant, vaccination is introduced to the susceptible individuals, susceptible individuals are recruited by birth or immigration, the treated individuals cannot transmit COVID-19 disease to the susceptible population, and by losing temporary immunity the recovered individuals become susceptible again. The compartments and parameters of the model are described in Table 1.

Qualitative analysis
In this section, we discuss the qualitative behavior of the model.

Positivity of the solution
Proof First, let us compute dS dt from model Eq. (2.1): Then applying the initial conditions and solving by using separation of variables, we get Next, from the second equation of (2.1) we get Finally, in a similar way, from the third, fourth, and fifth equations of model (2.1) we get Therefore our model is positive for all t ≥ 0.

Invariant region
To get the region at which the model is bounded, let us take the total population N(t) for all t, By differentiating Eq. (3.5) with respect to t we have In the absence of death due to COVID-19 (τ = 0 and ω = 0), Eq. (3.7) becomes Then solving Eq. (3.8) and applying the initial conditions, we get Thus our model is positively invariant in the region

Disease-free equilibrium point
In this section, we start with no COVID-19 present in the population of human, so that to obtain COVID-19-free equilibrium point, let us take I = T = R = 0. Now we have left with susceptible and vaccinated individuals from model Eq. (2.1): (3.11) Finally, we get the following points since I = 0, T = 0, and R = 0.

Basic reproduction number
In this section, we obtain both deterministic and stochastic reproduction numbers.

Basic reproduction number in the deterministic model
To get basic reproduction number, consider the newly infectious class of model Eq. (2.1): Let f be the rate of appearance of new infection, and let v be the rate of transfer of new infection into and out of the compartment. Then Then we obtain the Jacobian matrices of f and v (represented by F and V ) with respect to I: Then spectral radius (FV -1 ) of our model becomes Therefore our basic reproduction number in the deterministic model is . (3.14)
Now from Eq. (3.24) we observe that Therefore our disease-free equilibrium point is locally asymptotically stable if R D 0 < 1.

Local stability of disease-free equilibrium point in the stochastic model Theorem 3.3
If R S 0 < 1, then for any initial values of (S 0 , V 0 , Proof Let as take the Taylor series expansion of F(t, I(t)) = ln I(t). By Itô's formula Now integrating both sides of Eq. (3.25), we get where G(t) = t 0 β 3 dB 3 (t) is a martingale. By using the strong law of martingales we have lim t→∞ sup G(t) t = 0 almost surely. Then divide both sides of Eq. (3.27) by t: Now taking lim t→∞ sup of (3.28), we obtain Obviously, (μ + δ + τ ) > 0, and therefore Therefore our disease free-equilibrium point is locally asymptotically stable if and only if R S 0 < 1.

Global stability of disease-free equilibrium point
Theorem 3.4 If R D 0 < 1, then E 0 is globally asymptotically stable in .

Endemic equilibrium point
An endemic equilibrium point is a point in which the disease persists in the population, and we denote it as E 1 = S * , V * , I * , T * , R * different from zero. Now E 1 can be obtained by equating all system of model Eq. (2.1) to zero: Then after certain steps, we get , Therefore the endemic equilibrium point of equation system (2.1) is

Stability analysis of endemic equilibrium point
In this section, we study the stability of model (2.1) at the endemic equilibrium point. Proof Consider the Jacobian matrix at the endemic equilibrium point: The characteristic polynomial equation of J (E 1 ) is Hence by the Routh-Hurwitz stability criterion the endemic equilibrium point is locally asymptotically stable if This completes the proof.
Theorem 3.6 When R D 0 > 1, the endemic equilibrium E 1 is globally asymptotically stable.
Proof To show the global asymptotic stability at endemic equilibrium point, we define the following Lyapunov function: By calculating the derivatives of L(S * , V * , I * , T * , R * ) with respect to t we get Then simplifying Eq. (3.39), we have Therefore by Lasalle's invariant principle this indicates that the endemic equilibrium E 1 is globally asymptotically stable if U < V .

Sensitivity analysis and its interpretations
We determine the sensitivity of each parameter on the basic reproduction number of the model. We apply the following sensitivity index formula: where m i are the parameters of the basic reproduction number. Here The sensitivity analysis interpretation of our basic reproduction number is described as follows. The parameters that have negative sensitivity indices (δ, τ , θ ) have the effect of reducing the burden of COVID-19 from the community if the values of the two parameters are decreasing (which means that the basic reproduction number of the disease decreases as their parameter values decrease). Also, those parameters with positive sensitivity indices (σ 2 , ϕ, α) have an important role in the expansion of COVID-19 in the community if their values increase (this means that if their parameter values increase, then the secondary infection in the community increases). Therefore, so as to minimize the disease from the community, it is vital to decrease the positive indices and increase the negative indices. In the study of sensitivity, increasing the human mortality rate to control disease epidemic is not ethically acceptable, and hence we do not consider it. See Table 2 for more information.

Numerical simulation results and discussions
In this section, we obtain the numerical simulation results of the deterministic and stochastic models of COVID-19. To display the results, we used Maple 18 and MATLABsoftware. Now the graphs listed below demonstrate the results using the parameter values and initial conditions of the developed model. Here we have taken the parameter and variable values by assumption, estimation, and from recently published papers; they are listed in Table 3   α 0.0143 [19] ρ 0.0012 [23] μ 0.016 [19]

Deterministic and stochastic trends of the model
From the parameter values in Table 3 we obtained our basic reproduction numbers R D 0 = 0.0195522 and R S 0 = - 1.19667. We see that R D 0 < 1 and R S 0 < 1, which means that only susceptible population and vaccinated individuals are present and that infected, treated, and recovered populations are reduced to zero. This indicates that the model is asymptotically stable at R 0 < 1 for both deterministic and stochastic approaches, and this satisfies our theorem. This is verified numerically in Fig. 2.  Fig. 2 we observe that in both deterministic and stochastic approaches the number of infectious population will increase and susceptible population will decrease after a long period of time when the susceptible individuals contact with infected individuals. In addition, in the deterministic case the graph shows a smooth curve, and in the stochastic case, we see the zigzagging properties (which imply that there is a probability for susceptible individuals to be infected by COVID-19 when infected individuals contact with susceptible one). Therefore from this we can conclude that the stochastic approach is more advisable for such kinds of disease.

Effects of contact rate on infected individuals
As shown in Fig. 3, the simulation results of the effect of contact rate on infectious individuals I(t) are illustrated. This figure is displayed by varying the values of contact rate α and keeping the remaining parameter values unchanged. When the values of contact rate is large (α = 0.5), there is high possibility for the population to be infected by COVID-19, and for small values (α = 0.000001), there is low probability for the individual to be infected by COVID-19. Moreover, the amount of infectious individuals increases as the contact rate increases in both deterministic and stochastic approaches. Therefore we strongly advise for the concerned body that by decreasing contact with COVID-19 infected individuals it is possible to reduce the disease in the community.

Effects of treatment rate on infected individuals
We investigate the experimental results of impacts of treatment rate on infectious population I(t). As we see in Fig. 4, running numerical simulation results for deterministic approach is faster than for the stochastic approach, which clearly implies that the deterministic approach does not consider any probability (random properties) as the stochastic approach does. Besides, the infectious individuals become reduced by increasing treatment rate δ in both approaches. From this we can conclude that by treating the infectious individuals the infected population goes to treated one and COVID-19 will be eliminated from the community.

Effects of recovery rate on recovered individuals
In this subsection, we get the numerical simulation results of recovery rate δ on the amount of recovered individuals R(t). In Fig. 5 we display the simulation results by keeping different values of recovery rate (ρ = 0.00000035, ρ = 0.012, ρ = 0.03) and constant values of the other parameters. The figure clearly indicates that in the deterministic and stochastic approaches the increments of recovered individuals are obtained as the recovery rate in the individuals increase. Hence the greater the recovery rates, the more the individuals recover from COVID-19 in the community.

Discussions and conclusions
In this investigation, we presented both deterministic and stochastic models for a novel coronavirus COVID-19 disease dynamics. Since COVID-19 by nature has a rapid transmission and variation of spreading ratio in different environments, in this study the stochastic approach performed better and provided more accurate results. This study is devoted to implement coronavirus mathematical model by containing vaccination class. Vaccination is the best pharmaceutical controlling strategy for COVID-19 disease. The basic reproduction number is calculated for both deterministic and stochastic approaches by using the next generation matrix method. We analyzed the existence and stability of a disease-free equilibrium point. The disease-free equilibrium points are locally asymptotically stable when R D 0 < 1 and R S 0 < 1. Simulation results and analysis of the model are performed using combinations of all compartments by varying the contact and treatment rates of infected individuals and the recovery rate on recovered individuals.
From our numerical results we found that in both deterministic and stochastic approaches when susceptible population contacts with infected individuals, there is high probability to increase the number of infectious population and to decrease the number of uninfected (susceptible) population. Moreover, increasing the contact rate on population has an impact on the rate of spread of COVID-19 in the community. This means that when the values of contact rate α are large enough, there is high probability for the individuals to be infected by novel coronavirus COVID-19. In addition to this, infectious individuals are reduced by increasing the treatment rate δ. Hence vaccination strategies are regarded as the most effective measures to prevent and control rapid transmission of novel coronavirus COVID-19 in the community. Thus every citizen can take vaccine to minimize the spread of the disease. Therefore we strongly recommend other potential re-searchers to study different vaccination phases, like first and second doses of vaccinations (V 1 and V 2 ), and extensions by using optimal control strategies to modify our model.