A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate

: Infectious diseases remain a substantial public health concern as they are among the leading causes of death. Immunization by vaccination can reduce the infectious diseases-related risk of suffer-ing and death. Many countries have developed COVID-19 vaccines in the past two years to control the COVID-19 pandemic. Due to an urgent need for COVID-19 vaccines, the vaccine administration of COVID-19 is in the mode of emergency use authorization to facilitate the availability and use of vaccines. Therefore, the vaccine development time is extraordinarily short, but administering two doses is generally recommended within a speciﬁc time to achieve sufﬁcient protection. However, it may be essential to identify an appropriate interval between two vaccinations. We constructed a stochastic multi-strain SIR model for a two-dose vaccine administration to address this issue. We introduced randomness into this model mainly through the transmission rate parameters. We discussed the uniqueness of the positive solution to the model and presented the conditions for the extinction and persistence of disease. In addition, we explored the optimal cost to improve the epidemic based on two cost functions. The numerical simulations showed that the administration rate of both vaccine doses had a signiﬁcant effect on disease transmission.


Introduction
With a good understanding of the course of the epidemic, we will be in a better position to prevent or monitor communicable diseases, such as Coronavirus Disease 2019 (COVID- 19), and to help with the decision making for an effective control strategy. The mathematical epidemic models can forecast the epidemic trajectories and estimate the dynamics of disease transmission and recovery. This information can guide public health practitioners' responses, offer policymakers a reference for the best allocation of medical resources, and reduce the risk caused by the spread of the epidemic. Both deterministic and stochastic models have been proposed to model infectious diseases. Modern epidemic models are mostly constructed under the framework of the Susceptible-infected-recovered (SIR) model [1], and most of these models are presented in a deterministic system [2][3][4]. We refer you to [5] for a detailed presentation about the concepts and the related tools for deterministic models. Recently, many articles have introduced stochastic processes to epidemic models. For example, [6,7] brought in random disturbances to the transmission rate of the traditional Susceptible-infected-susceptible (SIS) model. Ref. [8] added random disturbances to the deterministic Susceptible-infectedrecovered-susceptible (SIRS) model and explored the impact of intervention strategies. Ref. [9] imported Levy noise in the SIS model. Quarantine and vaccination are the primary strategies against infectious diseases. Recent studies in the literature have considered models of quarantine [10][11][12][13][14] and have Table 1. Description of variables and parameters in model (1).

S(t)
the number of susceptible individuals at time t I k (t) the number of infectious individuals at time t who were infected by strain k R(t) the number of recovered individuals at time t V i (t) the number of individuals who received the ith dose of the vaccine at time t b the birth rate β k the transmission rate of strain k γ k the recovery rate of strain k δ the rate at which immunity disappears in recovered individuals q 0 the vaccination rate of the first dose among susceptible individuals q 1 the vaccination rate of the second dose among individuals who received first dose of vaccine ε 1 the rate at which immunity disappears among susceptible individuals receiving 1 dose of vaccine ε 2 the rate at which immunity disappears among susceptible individuals receiving 2 doses of vaccine µ k the death rate due to strain k µ N the natural death rate We consider a multi-strain SIR model as shown in (1). We assume that there are K types of strain and only one single vaccine available. The model includes susceptible population (S), the number of individuals who received only one dose of the vaccine (V 1 ), the number of individuals who received two doses of the vaccine (V 2 ), the number of people infected by the kth strain (I k ), and the number of individuals who were recovered (R). We assume that each susceptible person is infected with at most a single strain at a time. We present their relationship between the compartments in Figure 1, and introduce the following proposed model dI k (t) = [β k S(t)I k (t) − (µ N + µ k + γ k )I k (t)]dt, k = 1, 2, . . . , K,  We assume that the vaccine has a complete anti-epidemic effect on each strain, but its effect will disappear at a certain rate, returning the vaccinated person to the situation before the vaccination. Although this is different from reality, it may approximate the ineffectiveness of the vaccine against different strains.
Suppose that the number of total populations is ( ) = ( ) + ∑ =1 ( ) + We assume that the vaccine has a complete anti-epidemic effect on each strain, but its effect will disappear at a certain rate, returning the vaccinated person to the situation before the vaccination. Although this is different from reality, it may approximate the ineffectiveness of the vaccine against different strains.

The Equilibrium Point and the Basic Reproduction Number
We have some equilibrium points of our model on Λ, and its first point is the diseasefree equilibrium point P 0 = S 0 , I 0 1 , . . . , If the number of infected people is positive only for strain k and the rest is 0, the endemic equilibrium point is P * k = S * , I * k , R * , V * 1 , V * 2 , and k = 1, 2, . . . , K. Note that I * k = I * 1 , . . . , I * K means I * k > 0 and I * l = 0, ∀l = k. Using the next generation method, we can obtain the basic reproduction number. Set and Mathematics 2022, 10, 1804

of 22
The next generation matrix is Substituting S with S 0 , we have the basic reproduction number R 0k = β k S 0 µ N +µ k +γ k for strain k in the deterministic model, that is, k = 1, 2, . . . , K. This yields the basic reproduction number of model (1) with

Sensitivity Analysis of the Basic Reproduction Number to Parameters
Sensitivity analysis can be used to explore the effect on outcomes due to the changes in each parameter within a specified range in a mathematical model. In the studies of [26,27], the LHS/PRCC (Latin hypercube sampling/partial rank correlation coefficients) technique was used to analyze the sensitivity of parameters in R 0 . In this subsection, we use the same technique to find the parameters that have more influence on R 0 .
First, we set the probability functions of the parameters. We assign the birth rate b to have a uniform distribution between 0 and 5, and denote it as b~U(0, 5). Except for ε 2 , we set the rest of the parameters to have uniform distributions between 0 and 1. We assume to have a lower vaccine failure rate after two doses of the vaccine, and hence, ε 2 < ε 1 . We set ε 2~U (0, ε 1 ).
There are 3K + 7 parameters in model (1), and we set the number of simulations to 1000 times. In other words, we generate a 1000 × (3K + 7) LHS table as the input matrix and use R 0 as the output variable. Next, we calculate the PRCCs and test whether the coefficients are 0 s. We use the lhs package in R to generate the LHS table and use the sensitivity package to test the PRCCs. Note that the bootstrap method with repeated sampling 2000 times is used. The results with K = 3 are shown in Table 2 and Figure 2. Table 2. Partial rank correlation coefficients.
The results are significant at the 0.05 level (*) or the 0.01 level (**).
Except for ε 1 and q 1 , other parameters are significant, but these two parameters are directly related to q 0 and ε 2 . Both q 0 and q 1 are negatively correlated with R 0 , that is, the higher the vaccination rate, the smaller the R 0 . This is in line with the hypothesis that vaccines inhibit disease transmission. Conversely, ε 1 and ε 2 represent vaccine failure rates, which are positively correlated with R 0 . The most sensitive parameters are b and µ N . The ratio b/µ N represents the upper limit of the population. From Equation (5), this ratio has a definite linear relationship with the value of R 0 . set 2~U (0, 1 ).
There are 3K + 7 parameters in model (1), and we set the number of simulations to 1000 times. In other words, we generate a 1000 × (3K + 7) LHS table as the input matrix and use 0 as the output variable. Next, we calculate the PRCCs and test whether the coefficients are 0′s. We use the lhs package in R to generate the LHS table and use the sensitivity package to test the PRCCs. Note that the bootstrap method with repeated sampling 2000 times is used. The results with K = 3 are shown in Table 2 and Figure 2.   Except for 1 and 1 , other parameters are significant, but these two parameters are directly related to 0 and 2 . Both 0 and 1 are negatively correlated with 0 , that is, the higher the vaccination rate, the smaller the 0 . This is in line with the hypothesis that vaccines inhibit disease transmission. Conversely, 1 and 2 represent vaccine failure rates, which are positively correlated with 0 . The most sensitive parameters are b and . The ratio / represents the upper limit of the population. From Equation (5), this ratio has a definite linear relationship with the value of 0 .

The Stochastic Model
Following the spirit of the model in [6], we assume that the disease transmission rates fluctuate around some average values. Let Ω, F , {F t } t≥0 , P be a complete probability space, and β k ∼ β 0k + σ k dB k (t), where B k are independent standard Brownian motions, β 0k is the expected value of β k , and σ k are the amplitudes of β k , k = 1, 2, . . . , K. Thus, the proposed stochastic model becomes Note that model (12) is generalized from SIR or can be seen as a special case of the SVIRS model [28]. However, unlike the previous articles, we considered the impact of two vaccine doses on the epidemic under multi-strain infectious diseases.

The Asymptotic Behavior of the Stochastic Model
In this section, we study the mathematical properties of the model (12).

Existence of Unique Positive Solution
In Theorem 1, we show that there is a unique global and positive solution of model (12). That is, Λ is the positive invariant set of model (12) with probability 1.

Extinction
The threshold for the disappearance of the disease is an important property for the epidemic model. We define the basic reproductive number for each strain as The basic reproductive number for model (12) is To facilitate subsequent discussions, we define The following lemma is useful for our discussion. [29]). Let M = {M t } t≥0 be a real-valued continuous local martingale vanishing at t = 0 . Then

Lemma 1 (Strong Law of Large Numbers
and also Theorem 2. In model (12), consider the following two conditions: k . If one of the above conditions holds, then for an arbitrary initial vector θ(0) ∈ Λ, the solution In other words, the number of infectious individuals infected by the kth strain, I k , tends to zero exponentially with probability 1 for k = 1, 2, . . . , K.

Proof. See Appendix B.
Remark 1. The condition (b) of Theorem 2 implies R S 0k < 1. As we can express R S 0k as a function of S 0 , we have Note that f is the maximum of f . Then, That is, we can rewrite Theorem 2 as the following Theorem 3.

Persistence
In this subsection, we discuss the persistence of the disease. In this paper, we call a disease persistent if liminf

Theorem 4. For arbitrary initial vector
where k * = arg max where k * = arg min From Theorem 4, the value of ∑ K k=1 I k t should oscillate in the interval ξ, ξ if t is large enough. That is, some strains will persist, and after a certain period of time t, the average number of infected people will stay in this interval. We next discuss the conditions for the extinction or existence of the strains. Define Then, it is obvious that R S 0k * = R S 0 . We prove that, under specified conditions, strain k * is the uniquely persistent strain in Theorem 4.
Proof. See Appendix D.

Vaccination Rates
In this subsection, we discuss the utilities of the proposed model. This model treats vaccination as the primary disease preventative approach, and the relevant parameters in the model are the two-dose vaccination rates q 0 and q 1 . Using Equation (5) and letting where As people who have received their first dose of the vaccine are often more willing to take a second dose, we can simplify the assumptions of model (12). For example, we set q 0 = q 2 1 , then S 0 = A B +q 1 B +q 1 +q 2 1 +q 3 1 /C , and As discussed in Remark 1, R S 0k has a global maximum at S 0 = β 0k σ 2 k . It implies that if β k > σ 2 k S 0 , R S 0k is an increasing function of S 0 . Using (27), we have It indicates that if the left endpoint of Equation (29) is less than 1, we can eradicate diseases through increasing vaccination rates.
Let us consider a simple vaccine strategy: How to minimize the basic reproduction number with a finite cost? In order to encourage people to get vaccinated, the cost of incentive measures, in addition to general administrative work and vaccine costs, may also be added. We assume that the vaccine rate cost can be represented by an increasing and unbounded function. That is, where c 0 , c 1 > 0. Suppose the upper bound of the cost is C 0 . Then, the nonlinear programming problem is Another cost function to consider focuses on the number of infected individuals. We use the upper bound on the average number of infected individuals over a time interval (ξ) as an estimate of the number of infected individuals. Let Then, and The numerator of ∂ξ ∂S 0 is If , ξ is a strictly increasing function of S 0 . Then, ξ is a decreasing function of q 1 . We can consider a cost function as where c 2 > 0 can be seen as the cost of infected disease per person.

Numerical Examples
We illustrate the main results with some numerical examples. To simplify the discussion, we consider the scenario with only 3 strains, that is, we set K = 3. The first three examples are mainly to verify our theoretical results. Example 4 is to discuss the situations with different vaccination rates. In Examples 1 to 4, the initial values of the variables are given as S(0) = 6.75, R(0) = 0.1, V 1 (0) = 0.1, V 2 (0) = 0.05, and I k (0) = 1, k = 1, 2, 3. We provide Example 5 to further mimic the reality and present explanations in more detail.
We simulated the numbers of infected individuals for each of three strains, and we show the results in Figure 3. As the result of Theorem 2, the number of infected individuals approached 0 with time. We simulated the numbers of infected individuals for each of three strains, and we show the results in Figure 3. As the result of Theorem 2, the number of infected individuals approached 0 with time. We simulated the numbers of infected individuals for each of three strains. The results are shown in Figure 4, indicating that the number of infected individuals approached 0 with time. Then, we have max k∈{1,2,3} R S 0k = 0.2813 < 1 and max k∈{1,2,3} β 2 0k 2(µ N +µ k +γ k ) = 0.12 < 0.4 , k = 1, 2, 3. The model meets condition (b) of Theorem 2. As S 0 = 1.8750 and β 0k < S 0 σ 2 k , k = 1, 2, 3, it does not meet condition (a) of Theorem 2.
We simulated the numbers of infected individuals for each of three strains. The results are shown in Figure 4, indicating that the number of infected individuals approached 0 with time.    Then, the basic reproduction numbers are (R 01 , R 02 , R 03 ) = (3.1818, 2.5974, 1.2834), R S 01 , R S 02 , R S 03 = (3.1430, 2.5620, 1.2688), and R * 01 , R * 02 , R * 03 = (2.6610, 2.1212, 1.0873). We have k * = 1. As S 0 = 2.7273, it implies β 0k > S 0 σ 2 k , k = 1, 2, 3. We simulated the numbers of infected individuals of the three strains and we present the results in Figure 5.  We simulated the numbers of infected individuals of the three strains and we present the results in Figure 5.  (µ N +µ 1 +γ 1 ) . Obviously, only the first strain will persist, and the other two strains will go extinct.
Moreover, we assume cost model (31) with the costs c 0 = 1.2, c 1 = 1, and C 0 = 4. The optimal q 1 = 0.8851, and we have the optimal basic reproduction number R S 0 = 2.6804. Note that if no one is vaccinated at all, lim q 1 →0 R s 0 = 11.1458. If we assume cost model (36) with c 0 = 1.2, c 1 = 1, and c 2 = 5, then the optimal q 1 = 0.9044 and we have ξ = 3.8743. If no one is vaccinated at all, lim Example 4. We use the same setting for all parameters as that in Example 3 except for T and (q 0 , q 1 ). We set T = 100 and both q 0 and q 1 take 3 levels, 0.1, 0.5, and 0.9. The R S 0 's and ξ's related to different (q 0 , q 1 ) combinations are shown in Table 3. All R S 0 's and ξ's decrease as q 0 and q 1 increase. This echoes the sensitivity analysis in Section 2. It also means that the higher the vaccination rate, the higher the immune effect that can be achieved. We simulated the total number of infected individuals by three strains (∑ 3 k=1 I k (t)) and we present the results in Figure 6. The number of total infected individuals fluctuated within a fixed range after a certain time point. The level of its oscillation position is related to the vaccination rate.   Table 1 in [22].
As the 01 and 02 are both greater than 1, the number of infected patients gradually increases at the beginning and then decline after reaching a peak. The second strain slowly goes extinct, while the first strain continues oscillating between 25 and 50. 03 is less than 1, and it tends to go extinct from the beginning. This phenomenon can reflect that changes in the epidemic often return to a relatively flat fluctuation after a wave of peaks.

Conclusions
In this study, we construct a stochastic multi-strain SIR model based on the current phenomenon of COVID-19 and consider the situation of two-dose vaccine administration.  Table 1 in [22].
As the R S 01 and R S 02 are both greater than 1, the number of infected patients gradually increases at the beginning and then decline after reaching a peak. The second strain slowly goes extinct, while the first strain continues oscillating between 25 and 50. R S 03 is less than 1, and it tends to go extinct from the beginning. This phenomenon can reflect that changes in the epidemic often return to a relatively flat fluctuation after a wave of peaks.

Conclusions
In this study, we construct a stochastic multi-strain SIR model based on the current phenomenon of COVID-19 and consider the situation of two-dose vaccine administration. We use the basic reproduction number of this model, R S 0 , to establish the thresholds for disease extinction. In the discussions for persistence, we use another version of the basic reproduction number to find a lower bound on the number of infected patients. Although using the lower bound here may underestimate the number of infected patients, the result is mainly to guarantee the conditions for the persistent phenomenon. At the same time, we also point out the conditions under which some strains may go extinct in the competition of survival among multiple strains.
Finally, we introduce two cost models by leveraging the features of this model. That is, when minimizing the basic reproduction number or the number of infected individuals, we take the cost into consideration and provide the best two-dose vaccine administration rate. We also illustrate relevant theoretical results using numerical examples.
Our model has also a few limitations. The proposed model only considers the situation of one vaccine, and the failure mode of the vaccine is only expressed as the failure rate. Additionally, we do not discuss breakthrough infections and the potential for vaccines to reduce the disease severity and mortality. These are all directions that can be improved in the future.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Proof of Theorem 1. Let τ e be the explosion time. For any given initial vector θ(0) ∈ Λ, there is a unique solution θ(t) ∈ Λ for all t ∈ [0, τ e ). Let m 0 > 0 be a large value such that For any m > m 0 , we define the stopping time Obviously, τ m is increasing with m, and τ ∞ = lim m→∞ τ m ≤ τ e a.s. If we can prove τ ∞ = ∞ a.s., then τ e = ∞ a.s. and θ(t) ∈ Λ for all t ≥ 0 a.s.
If τ e = ∞ a.s. is invalid, there will exist a T > 0 and an > 0 such that It means there exists a large enough m 1 ≥ m 0 such that for all m ≥ m 1 , we have Define a function W : R K+4 By the Itô's formula, where Let Ψ m = {τ m ≤ T}, ∀m ≥ m 0 , and then P(Ψ m ) > . For all ω ∈ Ψ m , we have It implies which leads to the contradiction when m → ∞ , The proof is complete.

Appendix B
Proof of Theorem 2.