A Stochastic Analysis of a Siqr Epidemic Model With Short and Long-term Prophylaxis

. This paper aims to incorporate a high order stochastic perturbation into a SIQR epidemic model with transient prophylaxis and lasting prophylaxis. The existence and uniqueness of the global positive solution is proven and a stochastic condition in order to study the extinction of an infectious disease is established. The existence of a stationary distribution for the stochastic epidemic model is investigated as well. Numerical simulations are conducted to support our theoretical results and an example of application with COVID-19 data from Canada is used to estimate the transmission rate and basic reproduction number while constructing a model fitting the data.


Introduction
Investigating infectious disease dynamics has relied heavily on mathematical modeling.Since the seminal work of Kermack and McKendrick (1927), researchers employed deterministic and stochastic mathematical models to describe the unexpected behaviors and the evolution/recurrence of an epidemic.Due to its realistic approach, stochastic differential equations (SDEs) are excellent for describing such a phenomena.
In the literature, several stochastic epidemic models have been investigated such as Susceptible-Infected-Susceptible (SIS) (Wen et al., 2019), Susceptible-Infected-Removed (SIR) (Zhou et al., 2021a), Susceptible-Exposed-Infected-Removed (SEIR) (Jin and Lin, 2021), Susceptible-Infected-Removed-Susceptible (SIRS) (Lahrouz and Omari, 2013), and Susceptible-Infected-Removed-Infected (SIRI) (Zhou et al., 2022, El Fatini et al., 2018).Recent breakthroughs in SDEs allow us to include randomness into the description of biological events.Hence, numerous types of stochastic noise models have been introduced in epidemic modelling, including linear perturbations (Wen et al., 2019, Zhou et al., 2021a, Jin and Lin, 2021), Ornstein-Uhlenbeck processes (Wang et al., 2018, Cai et al., 2018) or incorporating jumps with Lévy noises (Berrhazi et al., 2017, Privault andWang, 2021).Recently, a new kind of stochastic perturbation has been developed, wherein time of a presence of a high volatility, the trending of a compartment realize some bursts describing an unknown behavior change.In fact, Liu and Jiang (2017), Liu et al. (2020b), Rajasekar et al. (2022) incorporated a non-linear perturbation for high environmental effects on the population dynamics and the infectious disease behavior.Liu and Jiang (2017) applied Lyapunov function technique on a stochastic SIR with non-linear perturbation to establish sufficient conditions for the existence of a unique ergodic stationary distribution of the model and derived sufficient conditions for an extinction of the disease.Also, Liu et al. (2020a) investigated an epidemic model with relapse and media coverage with higher order.Moreover, they derived sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of positive solutions to the dynamical system by establishing a suitable stochastic Lyapunov function.They also obtained adequate conditions for a wiping out of the infectious disease.In addition, to describe mathematically some different measures taken for managing an infectious disease epidemic, El Fatini et al. (2020) investigated a vaccination strategy in a stochastic epidemic models using a constant rate, where they established a stochastic threshold for the extinction and the existence of a stationary distribution for a persistent infectious disease with respect to delayed terms for the temporary immunity.Furthermore, several other works, e.g., Liu et al. (2018), Dai et al. (2022), Zhou et al. (2021b) studied other approaches to assess a vaccination strategy in a stochastic epidemic model, where they introduced a vaccinated compartment to explore the stochastic stability of the models.Moreover, some models, e.g., Liu et al. (2018) andEl Fatini et al. (2019), incorporated a loss of immunity by a perturbed transmission rate driven by a white noise.Finally, Caraballo et al. (2020), Zhang and Liu (2021), Zhang et al. (2017) included a quarantine strategy into stochastic epidemic models and they investigated their stochastic stability for the free disease and endemic cases.
The main purpose of this article is to show the existence and uniqueness of the stochastic Susceptible-Infected-Quarantined-Removed (SIQR) model with non-linear perturbation and short and long terms prophylaxis.Also, we improved results on stochastic R s 0 coefficients investigated by Liu and Jiang (2017) by including vaccination and non-pharmaceutical interventions (NPIs) and we support our theoretical results by numerical simulations and a case study using a real dataset.
The article is organized as follows: In Section 2, we describe the deterministic and stochastic epidemic models.In section 3, we present the main results for the stochastic model, i.e., the existence and uniqueness of the global positive solution (Section 3.1), and we find sufficient conditions for the extinction (Section 3.2) as well as for the existence of a stationary distribution for the persistence case of the infectious disease (Section 3.3).An example of application is presented in Section 4 using a case study of COVID-19 cases with the vaccination campaign in Canada, together with simulation results illustrating the long-term properties of our model using finite samples.Finally, in Section 5, we conclude with some research perspectives.The proof of all results are relegated to the Appendix.

Deterministic and stochastic models
Mathematical models such as El Fatini et al. (2020Fatini et al. ( , 2019) ) included a vaccination policy to a SIQR epidemic model to discuss the impact of quarantine and vaccination policy into a compartmental model, where S(t) describes the susceptible individuals, I(t) denotes the infected compartment, Q(t) denotes the quarantined individuals, and R(t) is the recovered compartment, including vaccinated individuals.
Our proposed model aims to add a higher order environmental perturbation with a general incidence function and a mitigation function describing NPIs.More precisely, the deterministic epidemic model is described by the following system of differential equations: where µ is the natural death rate in the population, i.e., for all compartments, β denotes the (unknown) transmission rate coefficient from susceptible to infected individuals, γ 1 denotes the vaccination rate of susceptible individuals, γ 2 is the recovery rate of the infected individuals, α 2 and α 3 represents the death rate for infected and quarantined individuals caused by infection complications, δ is the rate of infectious individuals who were quarantined, and γ 3 describes quarantined people who recovered from the infection.In case there is no quarantine, we set δ = Q 0 = 0.In addition, the general incidence function f (El Fatini et al., 2019, Lahrouz, 2015, Korobeinikov, 2007) x is non-increasing on [0, ∞), with f (0) = 0, and ϕ(0) = f ′ (0) = 1.Hence, ϕ(x) ≤ ϕ(0) = 1.These requirements are met by a number of well-known incidence rates, such as those indicating the effect of saturation (Capasso and Serio, 1978) on the infected individuals, where , or non-monotone incidence functions (Xiao and Ruan, 2007) to describe the psychological effects, where f (x) = I 1 + rx 2 , and in case of r = 0, we get the mass-action form with bilinear interactions f (x) = x.Since prophylaxis measures are not instantaneous and take time to be effective/implemented, the delayed effect must be taken into account.To this end, Betti et al. (2021) proposed to mitigate the effect of prophylaxis measures by using a function like M (t) = a + (1 − a)e −mt , where a ∈ (0, 1).Since M is always multiplied by β, we assume that M (t) ∈ (0, 1] for all t ≥ 0, with M (0) = 1 and lim t→∞ M (t) = a ∈ (0, 1).In what follows, this function is used to capture partially the effects of NPIs on the transmission rate, since it is continuous, with an initial point 1, and it converges to the parameter a, reflecting the value of an NPIs establishment by decision makers.
The asymptotic behavior of the deterministic system epidemic model ( 1) is determined by the re- More precisely, we have the following result, proven in Appendix , where Ī is the unique solution of , so there is a unique Next, to introduce randomness in the model, let W 1 , W 2 , W 3 , W 4 be independent Brownian motions starting at 0, defined on a complete probability space (Ω, F, P) and assume that σ ij represent volatility parameters, i ∈ {1, 2, 3, 4}, j ∈ {1, 2}, and constant rates µ i are replaced by the stochastic terms . More precisely, our proposed stochastic model is given by the following system of SDEs: A particular case of this stochastic model was investigated by Liu and Jiang (2017), where p = 0, M ≡ 1 and f (x) = x, without a quarantine compartment Q.

Main theoretical results
Throughout this section, we will consider the three compartments S(t), I(t) and Q(t) to study the behavior of the system of SDEs, since R(t) does not appear in the first three equations defining S, I, Q.Consequently, following (2), an investigation of compartments S, I, Q yields a solution for R.

Existence and uniqueness of the global positive solution.
In what follows, we first show that the solution of system (2) is global (i.e., without explosion in finite time) and non-negative.In Øksendal (2003), it is shown that a system of SDEs has a unique global solution for any given initial value, if the coefficients are Lipschitz and satisfy linear growth conditions.However, in practice, these conditions are too restrictive.Therefore, in (Mao, 2008, Theorem 3.5), the discontinuous functions are included in a stochastic differential equation by verifying the local Lipschitz condition and monotonicity conditions.Fulfilling the Has'minskii condition (Khasminskii, 2012) with a Lyapunov-like function, we get a sufficient condition for the existence and uniqueness.In fact, using the Lyapunov analysis method (Mao, 2008), it is simple to show that the solution of ( 2) is positive and global.The proof of the following result is given in Appendix B.2. Theorem 3.1.For any initial value (S 0 , I 0 , Q 0 ) ∈ R 3 + , there is a unique solution (S, I, Q) ∈ R 3 + of the stochastic system (2) on t ≥ 0 and the solution will remain in R 3 + with probability one.Remark 3.1.Note that our result holds also in the deterministic case (Ma et al., 2018) We recall from Miyahara (1972), the definition of a stochastic ultimately bounded process.
Definition 1.A process Z is said to be ultimately bounded if there exists a constant K such that for any (t, z) ∈ [0, ∞) × R d the following inequality holds: The proofs of the following theorem and its corollary are given in Appendices B.3-B.4.
The last inequality for E[T (t)] implies that the sequence T (n) is tight and as a result, there exists a subsequence n k so that T (n k ) converges in law.

Corollary 3.3. For any initial value
with probability one.

Extinction of the disease.
The extinction of infectious disease refers to the full and permanent lowering to zero of infected cases through deliberate efforts.In what follows, we investigate conditions to eliminate in a region or eradicate globally an infectious disease according to the stochastic system (2).To this end, set , and define , y > 0.
Further set f (y) = e g(y) (σ11+σ12y) 2 and π(y) = f (y)/C, so that π is a density.We are now in a position to state the extinction result, whose proof is given in Appendix B.5. Remark 3.2.Note that because of randomness, one can have extinction in the stochastic system while there is not extinction in the deterministic system, because one can have R 0 > 1 > Rs 0 , since 3.3.The existence of a stationary distribution.Persistent infectious diseases are identified as those in which the epidemic is not extinguished and remains in a population.For this matter, we investigate in this section sufficient conditions for the existence of a stationary distribution, under the assumption that M (t) = a for all t ≥ 0. To this end, for α and define Note that in particular, if . The proof of the following result is given in Appendix B.6.It extends and improve Theorem 2.1 in Liu and Jiang (2017), where, in their case M ≡ 1, f (x) = x, and their Řs Theorem 3.5.Assume that R s 0 (0) > 1.Then, for any initial value (S 0 , I 0 , Q 0 ) ∈ R 3 + , the system (2) has a unique stationary distribution π and it has the ergodic property.

Example of application and simulations
4.1.Case study of COVID-19 cases in Canada.In this case study, we consider data of COVID-19 disease spread in Canada from April 1, 2021 to May 30, 2021.We will try to predict the dynamics of COVID-10 for June 2021.Since there is no quarantine strategy, we set δ = 0 and we include a vaccinated population with two doses of Spikevax de Moderna, Vaxzevria d'AstraZeneca and Comirnaty de Pfizer-BioNTech or one dose of Janssen (Johnson & Johnson) COVID-19 vaccines to deduce the vaccination rate.The study period represents a wave that occurred in the first stages of the vaccination policy in Canada.We collected data from https://health-infobase.canada.ca,which provides precise information regarding the virus's spread throughout time and across the country.4.1.1.Parameter estimation.The interpretation of parameters is easier if the timescale is in years.In this case, if the observations are daily, then the difference h of time between observations should be h = 1 365 .Recall that for a deterministic or stochastic differential equation, the approximate solution We consider a data set of COVID-19 cases in Canada defined as Y T {y 0 , y 1 , . . ., y T }, where the observation is up to a finite horizon T and y t [S t , I t , R t ] ′ is a column vector in R 3×1 , which represents the daily observed values at time t for the susceptible (S), infected (I), and recovered individuals (R).We need to estimate the unknown vector of parameters θ in the epidemic model 1.To this end, we consider the predicted epidemic model ŷt (θ) defined by where ŷ0 (θ) = y 0 as an initial condition at time 0. Using this Euler method to approximate the solution, we calculate the quadratic cost Minimizing the quadratic cost J yields the non-linear least square estimator θ e , where In order to calculate the standard historical volatility where we look back over the historical data dynamics of the COVID-19 in Canada, we calculate u t defined as Doing so, we obtain the volatilities  1 shows the estimated parameters derived from fitting the studied models to the provided cumulative case data for Canada from April 1, 2021 to May 30, 2021.The predicted R = 0 1.0276 is, as expected, greater than 1, which means the disease will persist despite the vaccination deployment.In Figure 2, we used the estimated parameter β and the parameters deduced from the data to illustrate the patterns of the susceptible (S), infected (I), and recovered individuals (R).However, since this period represents the beginning of the vaccination policy the number of fully vaccinated population increased exponentially, which explains the alteration in the predicted pattern of susceptible and recovered compartments from the fitted data.This period illustrates high volatility since the SARS-CoV-2 B.1.617.2Delta variant was discovered in Canada in early April in British Columbia For this matter, we incorporate in Figure 3 a stochastic trajectory and higher order perturbation into the fitted deterministic trajectory from 60 days to predict the following 20 days including M (t) = a + (1 − a)e −mt , with a = 0.95 and m = 0.5 to describe NPIs.In Figure 4, we conducted 50 simulations to illustrate the probability density function of the infected population for the studied model.Therefore, the occurrence of infected cases in the population follows a stationary distribution and since R 0 > 1, the COVID-19 persists in the Canadian population regardless of the employment of the vaccine strategy.4.2.Vaccine deployment and quarantine policy.We recognize that there are many barriers to deploy a vaccination strategy because of human hesitancy and refusal towards different factors such as a mistrust of government and other institutions, conspiracy theories and quick process to test the vaccine.The major goal of the quarantine policy is to isolate the number of infected persons in the population in order to reduce interaction, resulting in a flattening of the spread curve, which leads to avoid a bigger pandemic wave.In Figure 5, we establish a quarantine strategy for 1 per cent of the daily infected population.A combination of a vaccination and quarantine policies with this rate will flatten the curve by reducing the probability of a contact between a susceptible and an infected individuals, and controlling the size of the infected population.The basic reproduction number value is 1.0478 with respect to a quarantine rate δ = 0.01, α 3 = α 2 and γ 3 = γ 2 .However, a significant proportion of the infected population shows no symptoms Long et al. (2020), and tracking accurately the infected individuals is not an easy task.Figure 6 shows a simulation the impact of doubling the vaccination rate to γ 1 = 0.0006, and relaxing the quarantine measure δ = 0.001.Therefore, the susceptible population decrease leading to a low probability of contact with infected individuals and it prevents also a possible

Discussion and concluding remarks
The objective of this research was to investigate a stochastic SIQR epidemic model driven by a non-linear perturbation with a short-term prophylaxis, such as social distancing and wearing masks, long-term prophylaxis like vaccination, and a general incidence function.The corresponding dynamical system was analyzed according to derived stochastic parameters to determine the extinction and the persistence of the infectious disease.This allowed us to explore Canadian data for COVID-19 cases from 1st April to 31th May in order to predict the trend of the pandemic in the following 20 days.Using least squares approach, we fitted the stochastic model and estimated the transmission rate and the basic reproduction numbers R 0 and R s 0 .We illustrated different scenarios to explore the utility of the combination of different strategies preventing and controlling the spread of the disease.Besides, Canada records above 50 percent of the maximum amount of moisture the air can hold, namely a high level of humidity percentage.For this matter, we explored a non-linear perturbation to describe the bursts in the infectious disease trending for the Canadian population in the spring season.Furthermore, it is important to note that the method proposed in this paper can be used to investigate other epidemic models, such as the SIHR model (Jiao and Huang, 2020).For future perspectives, it is intriguing to incorporate colored noise into the stochastic epidemic model (2), such as a continuous-time Markov chain.The motivation is that population dynamics may be affected by sudden environmental changes, such as temperature, humidity and so on.When switching between environments is quite often memoryless, the sudden-environmental changes can be modeled by a continuous-time Markov chain (Boukanjime et al., 2020).

Appendix A. Auxiliary results and notations
Suppose that Z = (Z 1 , . . ., Z d ) satisfies the following system of SDEs: where ã = σσ ⊤ , and Applying Itô's formula, and assuming z(t) ∈ S h , for all 0 ≤ t, we get ( 5) Next, let X(t) be a homogeneous Markov process in E d ⊂ R d described by the stochastic differential equation ( 6) Then there exists a stationary probability measure π for X so that Remark A.1.In Zhu and Yin (2007), to fulfill the Conditions (C1) and (C2), it suffices to show that there exists a bounded domain U with regular boundary and a non-negative C 2 -function V such that a(x) is uniformly elliptical in U and for any x ∈ E d \ U , LV (x) ≤ −c for some c > 0.
Taking τ k as in Theorem 3.1, letting k → ∞, and then multiplying both sides by e −µt , we get The latter is bounded in t since the right-hand side tends to A µ+ρ at t → ∞.Hence, according to the Definition 1, the process T (t) is stochastically ultimately bounded.

Figure 2 .
Figure 2. Data fitting of 60 days with Canadian COVID-19 cases and the prediction pattern for the next 20 days.pandemic wave for the Canadian population, which it can lead even to an extinction of the disease as R 0 = 0.5421 in absence of a new mutated variant of COVID-19 in the future.

Figure 3 .
Figure 3. Realistic , deterministic, stochastic and higher order stochastic trajectories to predict the next 20 days.
where a function b(t, z) = (b 1 (t, z), . . ., b d (t, z)) ⊤ defined for (t, z) ∈ [t 0 , ∞)×R d , and the matrix function σ(t, z) is a d × d matrix, b and σ are locally Lipschitz, and W = (W 1 , . . ., W d ) ⊤ is a d-dimensional Wiener process.Let Sh = {z ∈ R d : |z| < h}.The differential operator L associated with the law of the process Z, and acting on a function

Figure 4 .
Figure 4.The histogram of the probability density function of the infected population of the stochastic model without quarantine.
The following result corresponds to Theorem 4.1 inKhasminskii (2012).Lemma A.1.Suppose there exists a bounded open domain U ⊂ E d with regular boundary Γ having the following properties: (C1) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix a(x) is bounded away from zero.(C2) If x ∈ E d \ U , the mean time τ at which a path issuing from x reaches the set U is finite, and

Figure 5 .
Figure 5.A deployment of quarantine strategy in 1st April