Global asymptotic stability, extinction and ergodic stationary distribution in a stochastic model for dual variants of SARS-CoV-2.

Several mathematical models have been developed to investigate the dynamics SARS-CoV-2 and its different variants. Most of the multi-strain SARS-CoV-2 models do not capture an important and more realistic feature of such models known as randomness. As the dynamical behavior of most epidemics, especially SARS-CoV-2, is unarguably influenced by several random factors, it is appropriate to consider a stochastic vaccination co-infection model for two strains of SARS-CoV-2. In this work, a new stochastic model for two variants of SARS-CoV-2 is presented. The conditions of existence and the uniqueness of a unique global solution of the stochastic model are derived. Constructing an appropriate Lyapunov function, the conditions for the stochastic system to fluctuate around endemic equilibrium of the deterministic system are derived. Stationary distribution and ergodicity for the new co-infection model are also studied. Numerical simulations are carried out to validate theoretical results. It is observed that when the white noise intensities are larger than certain thresholds and the associated stochastic reproduction numbers are less than unity, both strains die out and go into extinction with unit probability. More-over, it is observed that, for weak white noise intensities, the solution of the stochastic system fluctuates around the endemic equilibrium (EE) of the deterministic model. Frequency distributions are also studied to show random fluctuations due to stochastic white noise intensities. The results presented herein also reveal the impact of vaccination in reducing the co-circulation of SARS-CoV-2 variants within a given population.

Introduction

The “severe acute respiratory syndrome coronavirus 2” (SARS-CoV-2), the cause of the Coronavirus Disease 2019 (COVID-19) was first reported in Wuhan, China. SARS-CoV-2 is a single-stranded-RNA virus with a genetic material which mutates in a very short time [24]. A large number of patterns have shown that only few mutations can lead to a more severe disease with higher transmissibility and infectivity [42]. Different variants of SARS-CoV-2 have been reported in the last two years, for example; the B.1.1.7 (alpha) variant detected in the UK, the “highly transmissible” and deadlier B.1.617.2 (delta) variant detected in India in December 2020 [24] and the recent B.1.1.529 (Omicron) variant first reported in South Africa. These variants have been categorized as “variants of concern (VOC)”, by the World Health Organization (WHO) and were detected in almost every country of the world [24]. The recently reported variants such as Delta and Omicron are linked with surge in COVID-19 cases where they are dominantly in circulation and “Detection and spread of the Omicron VOC in many more countries is expected. However, some medical experts have reported that the disease due to the Omicron VOC is less severe compared to other variants” [42]. Notwithstanding, out of several control measures against COVID-19, different vaccines have proven to be highly effective against severe illness with any of the variants. These include, but not limited to: the BNT162b2 Pfizer-BioNtech vaccine, the mRNA-1273 Moderna vaccine, Johnson and Johnson vaccine, and several others [16], [50]. These vaccines have high efficacy against different SARS-CoV-2 variants, including the Delta and Omicron variant [6]. We refer to [5], [35], [49] for epidemiological studies on the efficacy of available COVID-19 vaccines against VOCs.

Dual infection is a phenomenon where an individual is simultaneously infected with two or more strains of the same virus. It can affect host immune responses and results in an increase of the viral population. Reports strongly support the possibility of concurrent infection with dual variants of COVID-19 [1], [2], [47]. Scientists in Brazil identified two cases where people were simultaneously infected with two different variants of COVID-19 [1]. Both cases involved two young women who had typical mild-to-moderate flu-like symptoms and did not become severely ill or require hospitalization. In one case, the two variants identified had been circulating in Brazil since the beginning of the pandemic. In the other case, the person was simultaneously infected with both an older strain of the virus, and with the P.2 variant first identified in Rio de Janeiro. In addition, the case of an unvaccinated elderly woman in Belgium who was found to be infected with both the alpha and beta variants of COVID-19, has been confirmed [2]. She tested positive for the two variants on the same day and then developed rapidly worsening respiratory symptoms, which later led to her death. Furthermore, dual infections with both omicron and delta were found in immuno-competent and immuno-compromised patients living in different geographical areas [15], [46], [55]. In recent years, a number of cases when individuals were infected with more than one strain of HIV have been identified [20], [22], [53]. The findings on dual infections were reported for influenza viruses [33], the Epstein-Barr virus [56] and other viruses. Liu et al. [28] reported that a patient hospitalized in Iceland was infected by two SARS-CoV-2 subtypes simultaneously in early March 2020. One strain of the SARS-CoV-2 coronavirus was more aggressive, while the second strain was a mutation from the original version of the coronavirus that appeared in Wuhan, China. Hashim et al. [32] also reported double infection in all 19 analyzed samples, with most of the detected mutations genetically related. The authors discussed that co-infecting strains could compensate for the damaging effect of the truncated spike protein. However, the authors in [47] observed, when studying the different dynamics between two strains, that, one strain could replace the other after several days.

Most epidemic models are always influenced by environmental factors, such as precipitation, temperature, relative humidity, etc. [17], [36], [37], [41], [44]. For human disease associated epidemics, the nature of epidemic growth and spread is random due to the unpredictability in human-to-human contacts [48]. Thus, the variability and randomness of the environment is felt through the different states of the epidemic [51]. In epidemic dynamics, stochastic models may be more ideal in modeling epidemics in many circumstances [29]. For instance, stochastic models are able to incorporate randomness of infectious contacts occurring in the different infectious periods [13]. It has equally been shown that stochastic models can provide additional degree of realism as compared with their deterministic counterparts [54]. In particular, Allen et al. [8] showed that stochastic models best answer the question of disease extinction better than the deterministic counterpart. Herwaarden et al. [54] pointed out that an endemic equilibrium in a deterministic model could disappear in its corresponding stochastic system due to stochastic fluctuations. Furthermore, Nasell [34] revealed that stochastic models provide better approach to describe epidemics for a large range of realistic parameter values as compared with their deterministic equivalents.

Several mathematical models have investigated the dynamics of SARS-CoV-2 [7], [11], [25], [38], [39], [40] and its strains [10], [12], [21], [45]. Most of the multi-strain SARS-CoV-2 models do not capture randomness. As the dynamical behavior of most epidemics, especially SARS-CoV-2, is unarguably driven by random factors, it is thus very important to consider a stochastic vaccination co-infection model for two strains of SARS-CoV-2. In this paper, we highlight our contributions as follows:

We have developed a mathematical model for SARS-CoV-2 strains with vaccination, incorporating incident co-infection with both strains.

The stochastic model is analyzed for existence and uniqueness of solution.

The stochastic model is used to examine conditions for the extinction of both strains and their co-infection.

The conditions, for which the stochastic system fluctuates around the endemic equilibrium of the deterministic system are investigated. This was done with the help of appropriately defined Lyapunov functions.

The conditions under which a unique stationary distribution exist for the stochastic model are examined.

The deterministic and stochastic models are simulated to support the theoretical results.

The remaining part of the manuscript is organized as follows. The model which governs the dynamics of the disease is formulated in Section 2. Deterministic analyses are carried out in Section 3. The existence and uniqueness of non-negative solution in global sense, detailed analyses establishing the sufficient conditions for the extinction of both strains of the disease, asymptotic behavior of the stochastic model, as well as existence of a unique stationary distribution of the proposed model are all presented in Section 4. The numerical scheme and simulations to validate analytical results are shown in Section 5. At the end, in Section 6, the major outcomes of the study is summarized and further directions are given.

Preliminaries

In this section, we recall some basic concepts from stochastic calculus and some known theorems needed in the sequel. Throughout the paper, it is assumed that is a complete probability space with filtration . The following notations are introduced: Now, consider the d-dimensional stochastic differential equation [31]: with the initial condition where, stands for a d-dimensional standard Brownian motion defined on the given probability space. Let denote the family of all non-negative functions defined on which are continuously twice differentiable in . The differential operator of system (1) is defined by [31]: If acts on a function , then where, .

By Ito’s lemma [31], if , we have

Model formulation

At any time , the total population consists of the following epidemiological states: Susceptible individuals: , individuals infected with SARS-CoV-2 strain 1: , individuals infected with SARS-CoV-2 strain 2: , individuals co-infected with both strains: , individuals who have recovered from either or both strains of SARS-CoV-2: . From now onwards, strain 1 denotes the original strain of SARS-CoV-2, while strain 2 denotes other variants, such as delta and omicron variants, which are more transmissible than the original strain of SARS-CoV-2, and hence are the variants of concern (VoCs).

The model is formulated based on the following assumptions:

Individuals in the susceptible state can acquire SARS-CoV-2 strain 1 at the rate or strain 2 at the rate , where and stand for effective contact rates for incident infections.

As clinical reports affirm simultaneous infection with different variants of SARS-CoV-2 [1], [2], [28], [32], it is assumed that susceptible individuals can get co-infected with both strains from infected persons at the rate , where, is the effective contact rate for the transmission of dual strains of SARS-CoV-2.

Individuals already infected with SARS-CoV-2 strain 1 only can also acquire SARS-CoV-2 strain 2 at the rate , while those infected with SARS-CoV-2 strain 2 only can acquire strain 1 at the rate , with and denoting effective contact rates for infected individuals.

Susceptible individuals are vaccinated at the rate , and are assumed to have immunity against infection. However, due to the imperfect nature of the vaccine, this immunity is not life-long and can wane at the rate , upon which a vaccinated individual returns back to the susceptible class, where they can get infected with any of the SARS-CoV-2 variants.

Individuals in the susceptible group suffer natural death (just as those in other epidemiological states) at the rate .

Disease induced death is assumed for singly-infected and co-infected individuals at the rates and , respectively.

Parameters in the model are well described in Table 1 while the governing equations and flow charts describing transitions in the model are given in (2) and Fig. 1, respectively.

Table 1

Model parameters and variables.

Parameter	Description	Value	Source
Λ	Birth/recruitment rate	0.01138	[3]
γ	Natural death rate	176×365	[3]
ξ_1	Strain 1 contact rate	0.2944	[27]
ξ_2	Strain 2 contact rate	0.45	Assumed
ξ_12	Co-infection contact rate	0.42	Assumed
ψ_1	Strain 1 contact rate for infected persons	0.2944	[27]
ψ_2	Strain 2 contact rate for infected persons	0.45	Assumed
β_1	Strain 1 recovery rate	[130,13]/day	[14]
β_2	Strain 1 induced death rate	0.0214	[4]
α_1	Strain 2 recovery rate	[130,13]/day	[14]
α_2	Strain 2 induced death rate	0.0214	[4]
θ_1	Co-infection recovery rate	[130,13]/day	[14]
θ_2	Co-infection induced death rate	0.004	Assumed
η	Vaccination rate of SARS-CoV-2	0.01	Assumed
δ	Waning vaccine-induced immunity of SARS-CoV-2	0.01	Assumed

Fig. 1

The detailed flowcharts of SARS-CoV-2 transmission for systems (2), (3) respectively.

As, the dynamical behavior of most epidemics, such as SARS-CoV-2, is unarguably influenced by random factors, we have considered the corresponding stochastic two-strain co-infection epidemic model with vaccination described below: where, , and are the standard Gaussian white noise intensities, which represent the random/environmental factors that could influence the dynamics of the disease in each of the epidemiological states, while and denote independent standard Wiener processes, respectively.

Deterministic analysis

In this section, we analyze the deterministic system (2).

Deterministic basic reproduction number

The deterministic system’s disease free equilibrium (DFE) is given by The reproduction number of the model (2) is obtained by using the method in [52]. The transfer matrices are given by: The basic reproduction number of the model (2) is given by

, where , and are the associated reproduction numbers for SARS-CoV-2 strain 1, strain 2, and the co-infection of both strains, respectively and are given by:

Local asymptotic stability of the disease free equilibrium (DFE) of the model

The DFE, , of the model (2) is locally asymptotically stable (LAS) if , and unstable if .

The local stability of the model (2) is analyzed by the Jacobian matrix of the system (2) evaluated at the disease-free equilibrium, and is given by: where, yields the characteristic polynomial given by The eigenvalue are , and the solutions of the equations: Applying the Routh–Hurwitz criterion, the equations in (6) will all have negative real parts if and only if .

Endemic equilibrium points of the deterministic model

Boundary equilibria

When the reproduction number , the system (2) has three boundary endemic equilibria and as follows:

Strain 1 only:

Strain 2 only:

Co-infection only:

where:

Co-existence of endemic equilibria

In this section, we shall study the conditions for the co-existence of endemic equilibria of the deterministic system (2), and then will investigate the same numerically through simulations. We now prove the following result:

The deterministic system (2) has a family of co-existence endemic equilibria when .

Let, At endemic equilibrium (EE), the deterministic system (2) becomes, The steady state solutions is given by:

Adding all the classes at steady states (10), we have, where Substituting the steady state solutions into the force of infection , and simplifying, we obtain the following: Also, substituting the steady state solutions into the force of infection , and simplifying, we have: Finally, substituting the steady state solutions into the force of infection , and simplifying, we obtain the following: Since all the model parameters are positive, the model (2) exhibit co-existence endemic equilibria when the associated reproduction numbers: and . The epidemiological implication of this result is that, when the associated basic reproduction numbers of the deterministic system are greater than unity, then both strains and their co-infection become endemic in the population.

Stochastic analysis

The existence and uniqueness of solution, the global asymptotic behavior, and conditions for extinction and ergodic stationary distribution of the stochastic system are examined in this section.

Existence of a unique solution

In order to study the dynamical behavior, the first important question is whether the global solution exist or not. Moreover, for a model describing the dynamics of population, the nature of the value of the solution is also a matter of great interest. In this section, we show that the solution of stochastic system (3) is global and non-negative. It is known that, for a stochastic differential equation to have a unique global solution (that is, no explosion in a finite time) for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [9], [31].

Given any initial states , the random system (3) possesses a unique solution for , which will remain in with probability one, that is:

For any given initial conditions , there is a unique local solution for , with denotes the explosion time [31]. Suppose that is such that stay within . Then, each , consider the stopping time: where , where denotes the empty set. Note that, is increasing as . Thus, , so almost surely (a.s.). We need to show that a.s., so that and a.s for all .

Suppose that , then there exist a pair of constants and such that , so there exists an integer with Define function by Applying Ito’s lemma [31] on (16), we have On simplification, we have and hence we obtain that where Note that , and hence Thus, we have Integrating both sides of (18) from to , we obtain that If we take expectation on the both sides, we have Let . Then, by (15), we have . Note that, for every , we have at least or or or or or which is equivalent to or . Since So, Finally, we have where, is the indicator function of . If , then , becomes a contradiction. Thus, the only possibility is that , which completes the proof.

Extinction of the disease

In modeling the dynamics of any infectious disease, it is important to study conditions under which the disease will go into extinction or die out from the population. In this section, we will show that if the white noise is sufficiently large, then the solution of associated stochastic system (3) will become extinct with probability one.

To proceed, let us consider the following notation and results:

Strong Law of Large Numbers [31]

Let be continuous and real valued local martingale vanishing at and be its quadratic variation. Then

Let be a solution of (3) subject to , then a.s., Moreover, if , then

This Lemma’s proof is omitted as we can draw conclusion following the arguments similar to those given in Lemmas and in [58].

The threshold quantity for the stochastic system (3) can be written as where The following theorem gives the necessary conditions for the extinction of infection.

For any given initial value the solution

of system (3) has the following properties: If

, then SARS-CoV-2 strain 1 goes into extinction almost surely (a.s).

, then SARS-CoV-2 strain 2 goes to extinction a.s.

, then the co-infection of both strains goes to extinction a.s.

, then SARS-CoV-2 strain 1 dies out with probability 1.

, then SARS-CoV-2 strain 2 dies out with probability 1.

, then the co-infection of both strains die out with probability 1.

It means, if conditions (a) and (b) hold, then That is, the disease goes into extinction with probability .

Moreover,

Integrating the model (3), we have,

If formula is applied to second equation of (3), we have If we integrate Eq. (27) within , then we have This can be re-written as Division by gives that By the Strong Law of Large number, . As , by taking the limit superior on both sides of (28), we obtain This implies, . In a similar manner, it can be shown that which implies that . Note that, which gives .

Again, if Eq. (27) is integrated within and divided by , then we have Moreover, is continuous (locally) and . From Lemma 4.2 and , we have If , then Eq. (29) becomes Eq. (31) implies Likewise, if formula is applied to third equation of (3), then we obtain that Note that , is also a locally “continuous martingale” and . By Lemma 4.2 and , we have If , then Eq. (33) results in Eq. (35) gives If formula is applied to fourth equation of (3), then we obtain that Also, is a “locally continuous martingale” and . With Lemma 4.2 and , we have If is satisfied, then Eq. (37) becomes The Eq. (40) implies that Using 5th equation of (26), and by applying Eqs. (32), (36) and (40), as , we have From Eq. (26) and using Eqs. (32), (36) and (40), we obtain that Suppose the operator is given by Obviously as . So we can get Likewise, the last equation of (26) yields then we obtain from Eq. (44), if , From Eqs. (43) and (45), we obtain and

Stochastic asymptotic behavior of the perturbed model (3) around the endemic equilibrium (EE) of the deterministic system (2)

In this section, by constructing a suitable stochastic Lyapunov function, we study the asymptotic behavior of the stochastic system (3) around the endemic equilibrium (EE) of the deterministic system (3). Stochastic Lyapunov functions have been successfully constructed in several models to study the asymptotic behavior of the perturbed system around the endemic equilibrium of the deterministic systems [18], [30], [43]. In studying the dynamics of an epidemic system, when the disease will die out and when the disease will prevail in a population, are important questions to answer. If , then there is an endemic equilibrium for the deterministic system (2) but not for the stochastic system (3) as the stochastic system does not have the endemic equilibrium. Also, for , we have already shown the existence of endemic equilibria for the deterministic model (2). We shall now examine the conditions under which the solution of the perturbed system (3) fluctuates around the endemic equilibrium of the deterministic system (3) for a special case of the model. Due to strong non-linearity of the model, Lyapunov functions shall be constructed for a special case of the model when .

Let be the solution of (3) with initial condition . Also, let the reproduction number of the deterministic system (2) , (so that the endemic equilibrium exists). Then the solution of the stochastic system (3) has the asymptotic property: where,

For a special case of the model and when , the model has a co-existence endemic equilibria satisfying the following equations: Consider the stochastic Lyapunov function: where, Applying the Ito’s lemma [31], which can be written as: substituting the expressions at steady state (47), we obtain that which can also be written as, Applying the Young’s inequality [43]: for and , and using the fact , we have Also, applying the Ito’s Lemma to , we obtain that Using, and , we have Similarly, we have Adding Eqs. (49), (50), (51), (52) give On simplifying, we obtain that This can be simplified to obtain the following: with, On taking the expectation, we have Dividing by and taking limit as give The epidemiological implication of Theorem 4.3 is that the solution of the stochastic system (3) fluctuates around the positive endemic equilibrium of the deterministic system (2), as time as long as the white noise intensities are very small, since the difference between the solution of the stochastic system (3) and the endemic equilibrium of the deterministic system (2) is small to indicate that the disease will persist in the population.

Existence of ergodic stationary distribution

In epidemiological modeling, it is important to determine under what conditions the disease will persist in the population. For the deterministic model (2), it was shown that the boundary and coexistence endemic equilibria exist when the associated deterministic reproduction number is above unity. In Section 4.3, we also discussed the conditions under which the solution of the perturbed system (3) will fluctuate around the endemic equilibrium (EE) of the deterministic system (2). By applying results from Khas’minskii [26], it is now shown that there exist an ergodic stationary distribution which reveals that the disease will persist in the population for a long time. The aim of this section is to examine the existence of the solution to the perturbed model (3) which is a stationary Markov process. The stationary solution implies that the disease can be persistent and cannot die out in the population. First, we present some definitions and known results about the stationary distribution. Let represent a time-homogeneous regular Markov process in described by the following SDE: The corresponding diffusion matrix is defined as

The Markov process will have a unique ergodic stationary distribution whenever there exists a bounded open domain having boundary of regular type with the characteristics:

, matrix is strictly positive-definite.

, the mean time is finite where is the time required for moving from to , and for all compact subsets . Then, where represents a function integrable with respect to the measure .

Let us define the parameter Then the random model (3) possesses a unique stationary distribution having ergodicity whenever .

For , there exists a unique solution .

Moreover, the diffusion matrix of (3) is given as Suppose , then we have with, .

Thus, condition (a) of Lemma 4.3 holds.

We now need to define a function .

Let where and are positive constants which we shall determine. Using lemma, we have the following: Therefore, we have which implies that Let with Consequently, It implies that, Moreover, define where, , shall be determined later. Note that where . We now show that possesses the least value , which is unique.

If we consider partial derivatives of with respect to each variable, we have It can be observed that, at the function have the unique stagnation point . Furthermore, the Hessian matrix of the function at is given by Note that is positive-definite. Hence, has the least value.

Since is continuous, by Eq. (61), we have possesses a unique least value in .

Now, consider a function : Applying formula, we obtain that which leads to the following assertion where The next step is to define the set with for being negligible constants, to be determined later. For convenience, let us divide the domain into sub-domains as follows: Next, we show that in all twelve regions which ultimately implies that on .

Case. Suppose , then using Eq. (63), we have If we choose a sufficiently small so that the right hand side of the inequality in (64) is less than or equal to zero, then for .

Just as in the proof above, we conclude that for , , , , .

Case. Consider , then by Eq. (63), we get If we say, , for a very large positive value of and the least value of so that the right hand sides of the inequality in (65) is less than or equal to zero, then we have for very .

Just as in Case , we can obtain that for , ,

, , .

Hence, there exist such that Therefore, Let , and suppose that represents the time for a path to move from to Taking the integral on both sides of (66) within , taking expectation and then applying the formula by Dynkin’s [19], we get Since the function is not negative, therefore Following arguments similar to those given in the proof of Theorem 4.4, it can be easily seen that . Thus, the model (3) is regular. However, if we let and , then a.s. Hence, by using Fatou’s lemma [57], we have that Now, , where represents the subset of (which is compact); a direct fulfillment of condition b in Lemma 4.3. Thus, guaranteeing that the model (3) has unique stationary distribution.

Numerical scheme and simulations

The perturbed model (3) is simulated in this section. The numerical scheme used is based on the Milstein’s higher order method [23] and presented as follows: where, , represent the independent Gaussian random variables, with normal distribution , and being step size. denote the white noise values.

Numerical simulations are carried out in this section, to validate theoretical results discussed earlier. The initial conditions used for the simulations are assumed as follows: . Unless, otherwise stated, the values of parameters used are given in Table 1. In Fig. 2, simulations of the deterministic system (2) and perturbed system (3) are presented, when the contact rates and the stochastic white noise terms are respectively given by , , so that , , and the associated stochastic reproduction numbers are given by . The figure reveals that both strains and their co-infection go into extinction exponentially with unit probability. This also confirms the conclusions of Theorem 4.2. The biological implication is that this results in the elimination of both strains and their co-infection with unit probability. Both deterministic and stochastic models show agreement, as they converge to the DFE.

Fig. 2

Simulations of for the deterministic and stochastic models when the associated stochastic reproduction numbers are less than one. Parameter values used are: so that, .

From Theorem 4.4, we have established that, the random model (3) possesses a stationary distribution, which is unique. For the simulations using the contact rates and stochastic white noise terms: , , so that , it is observed that, for weak white noise intensities, the epidemic will remain and persist within the population. This is depicted in Fig. 3, where both strains and their co-infection will persist in average, thereby satisfying the results of Theorem 4.4. This also satisfies the conclusions of Theorem 4.3, that for small or negligible white noise disturbances, the solution of the perturbed system fluctuates around the endemic equilibrium of the deterministic system (2). Hence, the stochastic system (3) has a unique ergodic stationary distribution. The probability distribution histograms for the various epidemiological states of the model under this scenario are depicted in Fig. 4.

Fig. 3

Simulations of for the deterministic and stochastic models. Parameter values used are: , so that .

Fig. 4

The probability distribution histogram of and for the stochastic model (3), when the associated reproduction numbers are greater than one. Parameter values used are: , so that .

The simulations of the deterministic and perturbed systems for the case when the associated reproduction number is , are depicted in Fig. 5. The results reveal that both strains and their co-infection co-exist and persist within the population. This confirms the conclusions of Theorem 3.2, for the deterministic model (2). It is also observed that the random model (3) fluctuates around the co-existence endemic equilibrium when the associated stochastic model reproduction numbers for single infection and co-infection are greater than unity. Thus, in order to control the spread of different strains and their co-infection within the population, policies must ensure serious preventive efforts against the different variants.

Fig. 5

Simulations of for the deterministic and stochastic models, to numerically explore the coexistence of both strains and their co-infection. Parameter values used are: , so that, .

Impact of white noise and vaccination

Simulations of the perturbed model (3) to assess the impact of stochastic white noise intensities are presented in Figs. 6(a), 6(b) and 6(c). It can be observed that increasing the white noise intensities, hastens progression to extinction, for single and co-infected compartments. This simulation agrees with the conclusions of Theorem 4.3, confirming that for small white noise intensities, the solutions of the stochastic system fluctuates around the endemic equilibrium. However, for larger values of the white noise terms, the solution does not fluctuate around the EEP. This shows that sustained efforts in increasing stochastic disturbances through mass vaccination of susceptible individuals, adequate care and treatment for infected individuals could greatly reduce the spread and circulation of SARS-CoV-2 variants and their co-infection within the population. Simulations of the perturbed model (3) to assess the impact of vaccination are presented in Figs. 7(a), 7(b) and 7(c). It can be observed that increasing the vaccination rates have very high positive impact on the classes of individuals with single or dual strains of SARS-CoV-2. This shows that sustained efforts in vaccinating susceptible individuals could greatly reduce the spread and circulation of both variants and their co-infection within the population.

Fig. 6

Simulations of for the deterministic and stochastic models, to show the impact of stochastic white noise terms on the infected compartments, when , so that, .

Fig. 7

Simulations of for the deterministic and stochastic models, to show the impact of vaccination on the infected components, when , so that, .

Conclusion

In this work, we have presented a new stochastic model for two variants of SARS-CoV-2. The existence and the uniqueness of unique global solution of the stochastic model was shown. Using appropriately constructed Lyapunov functions, the conditions under which the solution of the perturbed system will fluctuate around the endemic equilibrium of the deterministic model were also derived. Stationary distribution and ergodicity for the new co-infection model were also established. Numerical simulations were carried out to validate theoretical results on extinction and persistence of SARS-CoV-2 variants within the population. We investigated the situations when both deterministic and stochastic associated reproduction numbers are below one and also when they are greater than one. Frequency distributions to show random fluctuations due to stochastic white noises were also presented. The simulation results also investigated the impact of vaccination on the dynamics of SARS-CoV-2 variants within the population.

Highlights of the simulations are as follows:

for the scenarios when the stochastic white noise terms, are larger than certain thresholds, that is: , , and the associated stochastic reproduction numbers are given by , both strains and their co-infection go into extinction exponentially with unit probability. This is presented in Fig. 2.

it was observed that, for weak white noise intensities, the solution of the stochastic system fluctuates around the endemic equilibrium (EE) of the deterministic model. More-over, it was shown that the random model (3) possesses a unique stationary distribution. For instance, when contact rates and stochastic white noise terms are: , , so that , both strains and their co-infection remain and persist within the population. This is depicted in Fig. 3, where both strains and their co-infection persist in average, thereby satisfying the results of Theorem 4.3, Theorem 4.4.

increasing the vaccination rates resulted in high positive population level impact on the classes of individuals with single or dual strains of SARS-CoV-2. Thus, enhancing efforts in vaccinating susceptible individuals could greatly reduce the spread and circulation of SARS-CoV-2 variants within the population. This is shown in Fig. 7

The stochastic model proposed in this research is based on the dynamics of different variants of SARS-CoV-2 with vaccination. The model did not capture cross-immunity between the variants. This could be an extension on the model, with more realistic assumptions on the emerging variants of concern. More-over, the emergence of different variants of SARS-CoV-2 warrants further studies on their co-infections with other diseases, such as Hepatitis B virus, tuberculosis, influenza, Malaria and other diseases. We could therefore, consider a robust stochastic model for variants of SARS-CoV-2 and co-infection with other diseases. This direction of research is still open to be explored. Also, due to insufficient information and reliable data about emerging variants of concern, we could not fit our model to real SARS-CoV-2 data. We hope to do this with more realistic and reliable information about the dynamics of different variants of concern.

CRediT authorship contribution statement

Andrew Omame: Conceptualization, Formal analysis, Methodology, Writing – original draft, Software. Mujahid Abbas: Conceptualization, Writing – review & editing, Supervision. Anwarud Din: Validation, Visualization, Formal analysis, Methodology, Writing – original draft, Writing – review & editing.