Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Factors Affecting the Use of Digital Mathematics Textbooks in Indonesia
Next Article in Special Issue
Robust Variable Selection for Single-Index Varying-Coefficient Model with Missing Data in Covariates
Previous Article in Journal
Organization Patterns of Complex River Networks in Chile: A Fractal Morphology
Previous Article in Special Issue
Graph Network Techniques to Model and Analyze Emergency Department Patient Flow
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Yen-Chang Chang
1,*,† and
Ching-Ti Liu
2,*,†
1
Center for General Education, National Tsing Hua University, Hsinchu City 300, Taiwan
2
Department of Biostatistics, School of Public Health, Boston University, Boston, MA 02118, USA
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2022, 10(11), 1804; https://doi.org/10.3390/math10111804
Submission received: 25 April 2022 / Revised: 17 May 2022 / Accepted: 23 May 2022 / Published: 25 May 2022

Abstract

:
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 suffering 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 specific time to achieve sufficient 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 significant effect on disease transmission.

1. 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-infected-recovered-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 incorporated the social network into models. However, due to the high cost of quarantine, vaccination is among the most effective and long-term approaches. Hence, there is considerable literature discussing infectious disease models with vaccines. For example, [15,16,17] introduced vaccination into SIS models. Ref. [18] discussed the vaccination under the SIR model from the perspective of optimal control. Ref. [19] presented a multi-strain infectious disease model allowing the virus to mutate into a new strain when there is vaccination in the SIR model.
The first reported cases of COVID-19 emerged in 2019 and has infected people around the world. The rapid spread of COVID-19 in the past two years has presented challenges to the international public health and economy. As of 17 April 2022, over 500 million confirmed cases and over six million deaths had been reported globally [20]. Many countries have developed COVID-19 vaccines to prevent the spread of the epidemic. There are several widely known and administered vaccines available on a large scale, including AstraZeneca (AZ), Moderna, Pfizer-BioNTech (BNT), Janssen (Johnson & Johnson), and Sputnik V. Over time, considerable discussions on vaccine-related topics have emerged. For example, [21] presented a systematic review. They provided an overview of clinical data from COVID-19 vaccine development and discussed the strategies to develop vaccines from different infectious diseases, with the aim to aid in developing therapeutic approaches for treatment and prevention of COVID-19. Ref. [22] proposed a new model to discuss the vaccine-resistant issue of COVID-19. Ref. [23] modified the SEIR compartmental model to forecast the possible development of the second COVID-19 wave in European countries. Ref. [24] proposed a multi-regional, hierarchical-tier model for understanding the complexity and heterogeneity of COVID-19 spread and control.
Due to an urgent need for COVID-19 vaccines, the vaccine administration of COVID-19 is in the mode of emergency use authorization requiring less data than typically needed in the early stage. The vaccine development time is therefore extraordinary short; however, to achieve sufficient protection, administering two doses are, in general, recommended within a specific time. Depending on which vaccine was administered, the recommended time-period between doses or the rate of administration may vary.
Therefore, it would be important to obtain an appropriate interval between two vaccinations, which can be then verified through clinical experiment. Recently, [25] proposed an extended SEAIRD compartmental model to describe the epidemic dynamics with both single-dose and double-dose vaccine administrations. They utilized the deterministic framework of nonlinear constrained optimal control problems to solve optimal vaccine administration problems. In addition, they only focused on the scenario with only one strain, which does not resemble reality. In this study, we hope to conduct a broader discussion on vaccination frequency through the construction of mathematical models. Additionally, due to a fast-spreading and higher mutation rate of the virus, there is already evidence to have more than one mutated strain of the new coronavirus. To take these into consideration, we will therefore construct a stochastic, multi-strain SIR model for a two-dose vaccine.
We structure the rest of this paper as follows. Section 2 presents the reasoning of our model and discusses its basic properties, such as equilibrium points and basic reproduction numbers of the deterministic model. Section 3 then focuses on its properties associated with the main stochastic model. Some numerical examples are shown to illustrate the results of our model in Section 4. We conclude this work in Section 5.

2. The Model

Table 1 lists the variables and parameters of our model.
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
{ d S ( t ) = [ b S ( t ) k = 1 K β k I k ( t ) ( μ N + q 0 ) S ( t ) + i = 1 2 ε i V i ( t ) + δ R ( t ) ] d t , d I k ( t ) = [ β k S ( t ) I k ( t ) ( μ N + μ k + γ k ) I k ( t ) ] d t ,   k = 1 , 2 , , K , d R ( t ) = [ k = 1 K γ k I k ( t ) ( μ N + δ ) R ( t ) ] d t , d V 1 ( t ) = [ q 0 S ( t ) ( μ N + ε 1 + q 1 ) V 1 ( t ) ] d t , d V 2 ( t ) = [ q 1 V 1 ( t ) ( μ N + ε 2 ) V 2 ( t ) ] d t .
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 N ( t ) = S ( t ) + k = 1 K I k ( t ) + i = 1 2 V i ( t ) + R ( t ) . It implies
d [ S ( t ) + k = 1 K I k ( t ) + R ( t ) + i = 1 2 V i ( t ) ] = [ b μ N N ( t ) k = 1 K μ k I k ( t ) ] d t [ b μ N N ( t ) ] d t
Then, we have N ( t ) N 1 , t 0 , where
N 1 = { b μ N , if   N ( 0 ) b μ N ; N ( 0 ) ,     if   N ( 0 ) > b μ N .
To simplify the following discussion, we assume N ( 0 ) b μ N . Next, we will discuss the stability of model (1) on an invariant set Λ :
Λ = { θ = ( θ 1 , , θ K + 4 ) + K + 4 | k = 1 K + 4 θ k b μ N } ,
where θ 1 = S , ( θ 2 , , θ K + 1 ) = I = ( I 1 , , I K ) , θ K + 2 = R , ( θ K + 3 , θ K + 4 ) = V = ( V 1 , V 2 ) .

2.1. The Equilibrium Point and the Basic Reproduction Number

We have some equilibrium points of our model on Λ , and its first point is the disease-free equilibrium point P 0 = ( S 0 , I 1 0 , , I k 0 , R 0 , V 1 0 , V 2 0 ) , where
{ S 0 = b μ N ( μ N + ε 2 ) ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( μ N + ε 1 + q 1 ) + q 0 ( μ N + ε 2 ) + q 0 q 1 , V 1 0 = b μ N q 0 ( μ N + ε 2 ) μ N 2 + μ N ( ε 2 + ε 1 + q 0 + q 1 ) + ( q 0 ε 2 + q 1 ε 2 + q 0 q 1 + ε 1 ε 2 ) , V 2 0 = b μ N q 0 q 1 μ N 2 + μ N ( ε 2 + ε 1 + q 0 + q 1 ) + ( q 0 ε 2 + q 1 ε 2 + q 0 q 1 + ε 1 ε 2 ) , I 1 0 = I 2 0 = = I K 0 = R 0 = 0 .
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
{ S * = μ N + μ k + γ k β k , V 1 * = q 0 ( μ N + μ k + γ k ) β 0 k ( μ N + ε 1 + q 1 ) , V 2 * = q 0 q 1 ( μ N + μ k + γ k ) β k ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) , R * = γ k I k * μ N + δ , I k * = b β k ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) μ N ( μ N + μ k + γ k ) [ μ N 2 + μ N ε 2 + q 0 q 1 ] β k ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( μ N + μ k + γ k γ k / ( μ N + δ ) )
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
F = diag ( β 1 S ,   β K S )
and
V = diag ( μ N + μ 1 + γ 1 ,   μ N + μ K + γ K )
The next generation matrix is
F V 1 = diag ( β k S μ N + μ k + γ k , k = 1 , 2 , , K )
Substituting S with S 0 , we have the basic reproduction number R 0 k = β k S 0 μ N + μ k + γ k for strain k in the deterministic model, that is,
R 0 k : = b μ N β k ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( μ N + μ k + γ k ) [ ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) + q 0 ( μ N + ε 2 ) + q 0 q 1 ] ,
k = 1 , 2 , , K . This yields the basic reproduction number of model (1) with
R 0 = max { R 0 k , k = 1 , 2 , , K }

2.2. 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.
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 .

2.3. The Stochastic Model

Following the spirit of the model in [6], we assume that the disease transmission rates fluctuate around some average values. Let ( Ω , , { t } t 0 , P ) be a complete probability space, and β k ~ β 0 k + σ k d B k ( t ) , where B k are independent standard Brownian motions, β 0 k is the expected value of β k , and σ k are the amplitudes of β k , k = 1 , 2 , , K . Thus, the proposed stochastic model becomes
{ d S ( t ) = [ b S k = 1 K β 0 k I k ( μ N + q 0 ) S + i = 1 2 ε i V i + δ R ] d t S k = 1 K σ k I k d B k ( t ) , d I k ( t ) = [ β 0 k S I k ( μ N + μ k + γ k ) I k ] d t + σ k S I k d B k ( t ) ,   k = 1 , 2 , , K , d R ( t ) = [ k = 1 K γ k I k ( μ N + δ ) R ] d t , d V 1 ( t ) = [ q 0 S ( μ N + ε 1 + q 1 ) V 1 ] d t , d V 2 ( t ) = [ q 1 V 1 ( μ N + ε 2 ) V 2 ] d t .
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.

3. The Asymptotic Behavior of the Stochastic Model

In this section, we study the mathematical properties of the model (12).

3.1. 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.
Theorem 1.
Given arbitrary initial vector θ ( 0 ) Λ , model (12) has a unique positive solution θ ( t ) Λ for all t 0 with probability 1. That is,
P ( θ ( t ) Λ ,   t 0 ) = 1 .  
Proof. 
See Appendix A. □

3.2. 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
R 0 k S : = R 0 k σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 ,
k = 1 , 2 , , K . The basic reproductive number for model (12) is
R 0 S : = max { R 0 k S ,   k = 1 , 2 , , K } .  
To facilitate subsequent discussions, we define
x t = x t = 1 t 0 t x ( s ) d s   and   x , x t = 1 t 0 t x 2 ( s ) d s .  
The following lemma is useful for our discussion.
Lemma 1
(Strong Law of Large Numbers [29]). Let M = { M t } t 0 be a real-valued continuous local martingale vanishing at t = 0 . Then
lim t M , M t = ,   a . s .     lim t M t M , M t = 0 ,   a . s .
and also
limsup t M , M t t < ,   a . s .     lim t M t t = 0 ,   a . s .
Theorem 2.
In model (12), consider the following two conditions:
(a)
R 0 k S = R 0 k σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 < 1  and  β 0 k > S 0 σ k 2 , or
(b)
β 0 k 2 2 ( μ N + μ k + γ k ) < σ k 2 .
If one of the above conditions holds, then for an arbitrary initial vector  θ ( 0 ) Λ , the solution θ ( t ) Λ  obeys
limsup t ln I k ( t ) t < 0 ,   a . s .
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 0 k S < 1 . As we can express  R 0 k S  as a function of S 0 , we have
f ( S 0 ) : = R 0 k S = β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 .
Note that  f ( β 0 k σ k 2 ) = β 0 k 2 2 ( μ N + μ k + γ k ) σ k 2  is the maximum of  f . Then, β 0 k 2 2 ( μ N + μ k + γ k ) < σ k 2  implies R 0 k S < 1 . That is, we can rewrite Theorem 2 as the following Theorem 3.
Theorem 3.
In model (12), consider the following two conditions:
(a)
R 0 k S = R 0 k σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 < 1  and  β 0 k > S 0 σ k 2 , or
(b)
β 0 k 2 2 ( μ N + μ k + γ k ) < σ k 2 ,
If one of the above conditions holds,then for arbitrary initial vector  θ ( 0 ) Λ , the solution θ ( t ) Λ  obeys  limsup t ln I k ( t ) t < 0 ,   a . s .
That is, I k tends to zero exponentially with probability 1 for  k = 1 , 2 , , K .

3.3. Persistence

In this subsection, we discuss the persistence of the disease. In this paper, we call a disease persistent if liminf t k = 1 K I k > 0 ,   a . s . holds. We define
R 0 k * : = R 0 k σ k 2 2 ( μ N + μ k + γ k ) ( b μ N ) 2 ,
k = 1 , 2 , , K . As R 0 k * R 0 k S , R 0 k * > 1 implies R 0 k S > 1 .
Theorem 4.
For arbitrary initial vector  θ ( 0 ) Λ ,
(a)
If max { R 0 k * } > 1 , we have
liminf t k = 1 K I k t ξ _ : = ( μ N + μ k ^ * + γ k ^ * ) ( R 0 k ^ * * 1 ) β 0 k ^ * S 0 b ( ( μ N + μ k ˜ * + γ k ˜ * ) δ γ k ˜ * μ N + δ )   a . s . ,
where  k ˜ * = arg   max k : R 0 k * > 1 { ( μ N + μ k + γ k ) δ γ k μ N + δ }  and  k ^ * = arg   max k : R 0 k * > 1 { ( μ N + μ k + γ k ) ( R 0 k * 1 ) β 0 k } .
(b)
If there are some persistent strains with  β 0 k 0 > S 0 σ k 0 2  and  R 0 k 0 * > 1 , then
limsup t k = 1 K I k t ξ ¯ : = ( μ N + μ k ^ 0 + γ k ^ 0 ) ( R 0 k ^ 0 S 1 ) ( β 0 k ^ 0 S 0 σ k ^ 0 2 ) S 0 b ( ( μ N + μ k ˜ * + γ k ˜ * ) δ γ k ˜ * μ N + δ )   a . s . ,
where k ˜ * = arg   min k : R 0 k * > 1 { ( μ N + μ k + γ k ) δ γ k μ N + δ } and k ^ 0 = a r g   m a x k 0 : R 0 k 0 S > 1 { ( μ N + μ k 0 + γ k 0 ) ( R 0 k 0 S 1 ) ( β 0 k 0 S 0 σ k 0 2 ) } .
Proof. 
See Appendix C. □
From Theorem 4, the value of k = 1 K 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
k ¯ * = arg max k : R 0 k S > 1 R 0 k S .
Then, it is obvious that R 0 k ¯ * S = R 0 S . We prove that, under specified conditions, strain k ¯ * is the uniquely persistent strain in Theorem 4.
Theorem 5.
Suppose the conditions of Theorem 4 all hold. Assume that
lim t ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t = 0 ,   a . s .
and there are some  k k ¯ *  with  R 0 k ¯ * S > R 0 k S > 1 . If
(a)
β 0 k ¯ * σ k ¯ * 2 β 0 k σ k 2  and  σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) > σ k 2 ( μ N + μ k + γ k ) , or
(b)
β 0 k ¯ * σ k ¯ * 2 < β 0 k σ k 2  and  σ k 2 ( μ N + μ k + γ k ) > ( R 0 k S 1 R 0 k ¯ * S 1 ) σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) ,
we have  limsup t ln I k ( t ) ln I k ( 0 ) t < 0 ,   a . s . , i.e.,  lim t I k ( t ) = 0 ,   a . s .
Proof. 
See Appendix D. □

3.4. 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 A = b μ N ,   B = μ N 2 + μ N ( ε 2 + ε 1 ) + ε 1 ε 2 , and C = ε 2 + μ N , then
S 0 = A B + q 1 ( ε 2 + μ N ) B + q 1 ( ε 2 + μ N ) + q 0 ( ε 2 + μ N ) + q 0 q 1 = A B + q 1 B + q 1 + q 0 + q 1 q 0 / C ,
where B = B C = μ N 2 + μ N ( ε 2 + ε 1 ) + ε 1 ε 2 ε 2 + μ N = μ N + ε 1 . 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 1 2 , then S 0 = A B + q 1 B + q 1 + q 1 2 + q 1 3 / C , and
S 0 q 1 = A ( 1 B + q 1 + q 1 2 + a q 1 3 C ( B + q 1 ) ( 1 + 2 q 1 + 3 q 1 2 C ) ( B + q 1 + q 1 2 + q 1 3 C ) 2 ) = A ( 2 B q 1 + ( 3 B C + 1 ) q 1 2 + 2 q 1 3 C ( B + q 1 + q 1 2 + q 1 3 C ) 2 ) < 0 .
As discussed in Remark 1, R 0 k S has a global maximum at S 0 = β 0 k σ k 2 . It implies that if β k > σ k 2 S 0 , R 0 k S is an increasing function of S 0 . Using (27), we have
R 0 k S q 1 = R 0 k S S 0 S 0 q 1 = A μ N + μ k + γ k β k σ k 2 S 0 2 B q 1 + 3 B C + 1 q 1 2 + 2 q 1 3 C B + q 1 + q 1 2 + q 1 3 / C 2 < 0 .
R 0 k S is a decreasing function of q 0 . For q 1 ( 0 , 1 ) , we have S 0 ( b μ N μ N + ε 1 + 1 μ N + ε 1 + 2 + 1 / ( ε 2 + μ N ) , b μ N ) and
R 0 k S A μ N + μ k + γ k × β 0 k B + 1 B + 2 + 1 C σ k 2 2 A B + 1 B + 2 + 1 C 2 , β 0 k σ k 2 2 A , k = 1 , 2 , , K .
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,
C 0 ( q 1 ) = c 0 ln ( 1 q 0 ) c 1 ln ( 1 q 1 ) = c 0 ln ( 1 q 1 2 ) c 1 ln ( 1 q 1 ) ,
where c 0 ,   c 1 > 0 . Suppose the upper bound of the cost is C 0 . Then, the nonlinear programming problem is
minimize     R 0 S subject   to c 0 ln ( 1 q 1 2 ) c 1 ln ( 1 q 1 )   C 0   and   0 q 1 1
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
D = 1 b ( ( μ N + μ k ˜ * + γ k ˜ * ) δ γ k ˜ * μ N + δ ) .
Then,
ξ ¯ = ( μ N + μ k ^ 0 + γ k ^ 0 ) ( R 0 k ^ 0 S 1 ) ( β 0 k ^ 0 S 0 σ k ^ 0 2 ) S 0 b ( ( μ N + μ k ˜ * + γ k ˜ * ) δ γ k ˜ * μ N + δ ) = β 0 k S 0 σ k 2 2 ( S 0 ) 2 ( μ N + μ k ^ 0 + γ k ^ 0 ) ( S 0 β 0 k ^ 0 S 02 σ k ^ 0 2 ) D
and
ξ ¯ S 0 = 1 2 β 0 k S 02 σ k ^ 0 2 + ( μ N + μ k ^ 0 + γ k ^ 0 ) ( β 0 k ^ 0 2 S 0 σ k ^ 0 2 ) ( S 0 β 0 k ^ 0 S 02 σ k ^ 0 2 ) 2 D .
The numerator of ξ ¯ S 0 is
1 2 β 0 k σ k ^ 0 2 S 02 2 D σ k ^ 0 2 S 0 + ( μ N + μ k ^ 0 + γ k ^ 0 ) β 0 k ^ 0 = 1 2 β 0 k σ k ^ 0 2 ( S 0 2 ( μ N + μ k ^ 0 + γ k ^ 0 ) β 0 k ) 2 2 D ( μ N + μ k ^ 0 + γ k ^ 0 ) 2 σ k ^ 0 2 + ( μ N + μ k ^ 0 + γ k ^ 0 ) β 0 k ^ 0 .
If β 0 k ^ 0 > 2 ( μ N + μ k ^ 0 + γ k ^ 0 ) σ k ^ 0 2 , ξ ¯ is a strictly increasing function of S 0 . Then, ξ ¯ is a decreasing function of q 1 . We can consider a cost function as
C 1 ξ ¯ , q 1 = c 2 ξ ¯ c 0 ln 1 q 1 2 c 1 ln 1 q 1 .
where c 2 > 0 can be seen as the cost of infected disease per person.

4. 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.
Example 1.
We set the study duration T = 10 and the gap of simulation Δ t = 0.01 . Assume the parameters are given by b = 1 ,   σ k 2 = 0.001 ,   k = 1 , 2 , 3 , ( β 01 , β 02 , β 03 ) = ( 0.1 , 0.08 , 0.05 ) , ( μ 1 , μ 2 , μ 3 ) = ( 0.2 , 0.1 , 0.05 ) , μ N = 0.1 , ( q 0 , q 1 ) = ( 0.8 , 0.6 ) , ( ε 1 , ε 2 ) = ( 0.3 , 0.1 ) , and δ = 0.1 .
Then, we have m a x k { 1 , 2 , 3 } R 0 k S = 0.2941 < 1 and S 0 = 2.3810 . It implies β 0 k > S 0 σ k 2 and β 0 k 2 2 ( μ N + μ k + γ k ) > 0.001 , k = 1 , 2 , 3 . The model meets condition (a) of Theorem 2, but not condition (b) of Theorem 2.
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.
Example 2.
We set T = 10 and Δ t = 0.01 . Assume the parameters are given by b = 1 ,   σ k 2 = 0.4 ,   k = 1 , 2 , 3 , ( β 01 , β 02 , β 03 ) = ( 0.6 , 0.5 , 0.4 ) , ( μ 1 , μ 2 , μ 3 ) = ( 0.4 , 0.4 , 0.4 ) , μ N = 0.1 , ( q 0 , q 1 ) = ( 0.8 , 0.6 ) , ( ε 1 , ε 2 ) = ( 0.9 , 0.5 ) , and δ = 0.1 .
Then, we have m a x k { 1 , 2 , 3 } R 0 k S = 0.2813 < 1 and m a x k { 1 , 2 , 3 } β 0 k 2 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 β 0 k < S 0 σ k 2 , 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.
Example 3.
We set T = 300 and Δ t = 0.01 . Assume the parameters are given by b = 1 ,   ( σ 1 2 , σ 2 2 , σ 3 2 )   = ( 0.0025 , 0.0020 , 0.0020 ) , ( β 01 , β 02 , β 03 ) = ( 0.28 , 0.20 , 0.24 ) , ( μ 1 , μ 2 , μ 3 ) = ( 0.04 , 0.01 , 0.01 ) , ( γ 1 , γ 2 , γ 3 ) = ( 0.1 , 0.1 , 0.4 ) ,   μ N = 0.1 , ( q 0 , q 1 ) = ( 0.64 , 0.8 ) ,  ( ε 1 , ε 2 ) = ( 0.3 , 0.1 ) , and δ = 0.1 .
Then, the basic reproduction numbers are ( R 01 , R 02 , R 03 ) = ( 3.1818 , 2.5974 , 1.2834 ) , ( R 01 S , R 02 S , R 03 S ) = ( 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 β 0 k > S 0 σ k 2 , k = 1 , 2 , 3 .
We simulated the numbers of infected individuals of the three strains and we present the results in Figure 5.
Note that σ 2 2 β 01 σ 1 2 β 02 = 1.1200 > 1 , σ 1 2 ( μ N + μ 1 + γ 1 ) = 0.0104 > 0.0095 = σ 2 2 ( μ N + μ 2 + γ 2 ) , σ 3 2 β 01 σ 1 2 β 03 = 0.9333 < 1 , and σ 3 2 ( μ N + μ 3 + γ 3 ) = 0.0039 > 0.0013 = ( R 03 S 1 R 01 S 1 ) σ 1 2 ( μ 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 0 S = 2.6804 . Note that if no one is vaccinated at all, lim q 1 0 R 0 s = 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 q 1 0 ξ ¯ = 5.9681 .
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 0 S ’s and ξ ¯ ’srelated to different ( q 0 , q 1 ) combinations are shown inTable 3. All R 0 S ’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 ( k = 1 3 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.
Example 5.
In order to mimic the actual COVID-19 infection situation, we consult the parameter information in Table 1 in [22] for simulation. First, we set b = 10 and μ_N = 0.01 to meet the value of N * in the table, and we use the data in the table to set δ = 1 180 , μ k = 0.01 × 1 14 , and γ k = 0.99 × 1 14 , k = 1 , 2 , 3 . We assume that people who received one dose of the vaccine had twice the rate of loss of immunity compared with those who received two doses. Then, we set ε 2 = 1 180 and ε 1 = 1 90 . As the vaccination rate is between 0.001 and 0.015, we set ( q 0 , q 1 ) = ( 0.01 ,   0.015 ) . As the transmission rate is between 0.055 N and 0.18 N in the table, we assign ( β 1 , β 2 , β 3 ) = 1 1000 × ( 0.18 ,   0.15 , 0.06 ) and ( σ 1 , σ 2 , σ 3 ) = 1 1000 × ( 0.12 ,   0.10 , 0.04 ) .
Then, the basic reproduction numbers are ( R 01 , R 02 , R 03 ) = ( 1.4317 , 1.1931 , 0.4772 ) , ( R 01 S , R 02 S , R 03 S ) = ( 1.3946 , 1.1673 , 0.4731 ) , and ( R 01 * , R 02 * , R 03 * ) = ( 1.3433 , 1.1317 , 0.4674 ) .
We simulated the numbers of infected individuals of the three strains with T = 100 , and the initial conditions are given as S ( 0 ) = 967 , R ( 0 ) = V 1 ( 0 ) = V 2 ( 0 ) = 1 , and I k ( 0 ) = 10 , k = 1 , 2 , 3 . The results are presented in Figure 7.
As the R 01 S and R 02 S 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 03 S 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.

5. 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 0 S , 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.

Author Contributions

Conceptualization, Y.-C.C. and C.-T.L.; methodology, Y.-C.C. and C.-T.L.; validation, Y.-C.C. and C.-T.L.; formal analysis, Y.-C.C. and C.-T.L.; investigation, Y.-C.C. and C.-T.L.; writing—original draft preparation, Y.-C.C.; writing—review and editing, C.-T.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially supported by the Ministry of Science and Technology of Taiwan under the grant MOST 109-2118-M007-002.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

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 θ ( 0 ) A m 0 , where
A m 0 = { θ ( t ) Λ : 1 m 0 < min 1 i 2 K + 4 θ i ( t ) ,   max 1 i 2 K + 4 θ i ( t ) < m 0 } .
For any m > m 0 , we define the stopping time
τ m = inf { t [ 0 , τ e ) :     θ ( t ) A m } .
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
P ( τ e T ) > ϵ .
It means there exists a large enough m 1 m 0 such that for all m m 1 , we have
P ( τ m T ) > ϵ .
Define a function W : + K + 4 + { 0 } by
W ( θ ) = ( S 1 ln S ) + k = 1 K ( I k 1 ln I k ) + ( R 1 ln R ) + i = 1 2 ( V i 1 ln V i ) .
By the It o ^ ’s formula,
d W = ( 1 1 S ) d S + 1 2 S 2 ( d S ) 2 + k = 1 K [ ( 1 1 I k ) d I k + 1 2 I k 2 ( d I k ) 2 ] + ( 1 1 R ) d R + i = 1 2 [ ( 1 1 V i ) d V i + 1 2 V i 2 ( d V i ) 2 ] = L W d t + ( 1 1 S ) [ S k = 1 K σ k I k d B k ] + k = 1 K ( 1 1 I k ) [ S σ k I k d B k ] = L W d t + ( k = 1 K ( I k S ) σ k d B k ) ,
where
L W = ( 1 1 S ) [ b S k = 1 K β 0 k I k ( μ N + q 0 ) S + i = 1 2 ε i V i + δ R ] + 1 2 S 2 ( S 2 k = 1 K ( σ k I k ) 2 )                     + k = 1 K [ ( 1 1 I k ) ( S β 0 k I k ( μ N + μ k + γ k ) I k ) ] k = 1 K 1 2 I k 2 S 2 ( σ k I k ) 2                     + ( 1 1 R ) [ k = 1 K γ k I k ( μ N + δ ) R ] + ( 1 1 V 1 ) [ q 0 S ( μ N + ε 1 + q 1 ) V 1 ] + ( 1 1 V 2 ) [ q 1 V 1 ( μ N + ε 2 ) V 2 ]             [ b + i = 1 2 ε i V i + δ R + k = 1 K β 0 k I k + ( μ N + q 0 ) ] + 1 2 k = 1 K ( σ k 2 I k ) 2                     + [ k = 1 K γ k I k + ( μ N + δ ) ] + k = 1 K [ S β 0 k I k + ( μ N + μ k + γ k ) ]                           + k = 1 K 1 2 S 2 σ k 2 + [ q 0 S + ( μ N + ε 1 + q 1 ) ] + [ q 1 V 1 + ( μ N + ε 2 ) ]               b + ( K + 4 ) μ N + δ + q 0 + q 1 + ε 1 + ε 2 + k = 1 K ( μ k + γ k )                     + b μ N ( i = 1 2 ε i + δ + k = 1 K β 0 k + k = 1 K γ k + q 0 + q 1 ) + ( b μ N ) 2 k = 1 K ( β 0 k + σ k β 2 ) M ,
for all 0 < t < T . Then, we obtain
E W ( θ ( T τ m ) ) W ( θ ( 0 ) ) + M T < .
Let Ψ m = { τ m T } , m m 0 , and then P ( Ψ m ) > ϵ . For all ω Ψ m , we have
W ( θ ( τ δ , ω ) ) ( m 1 ln m ) ( 1 m 1 + ln m ) .
It implies
W ( θ ( 0 ) ) + M T E W ( θ ( T τ m ) χ ( ω Ψ δ ) ) ϵ [ ( m 1 ln m ) ( 1 m 1 + ln m ) ] ,
which leads to the contradiction when m ,
> W ( θ ( 0 ) ) + M T .
The proof is complete. □

Appendix B

Proof of Theorem 2.
(a)
Using model (12) and Equation (16), we have
{ S ( t ) S ( 0 ) t = b k = 1 K β 0 k S I k t ( μ N + q 0 ) S t + i = 1 2 ε i V i t + δ R t                                                                                                                                                                           1 t k = 1 K σ k 0 t S I k d B k ( s ) , I k ( t ) I k ( 0 ) t = β 0 k S I k t ( μ N + μ k + γ k ) I k t + 1 t σ k 0 t S I k d B k ( s ) ,   k = 1 , 2 , , K , R ( t ) R ( 0 ) t = k = 1 K γ k I k t ( μ N + δ ) R t , V 1 ( t ) V 1 ( 0 ) t = q 0 S t ( μ N + ε 1 + q 1 ) V 1 t , V 2 ( t ) V 2 ( 0 ) t = q 1 V 1 t ( μ N + ε 2 ) V 2 t .
According to system (A10) of equations, we have
( S ( t ) S ( 0 ) t + k = 1 K I k ( t ) I k ( 0 ) t ) + ε 1 ( μ N + ε 2 ) + ε 2 q 1 ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) V 1 ( t ) V 1 ( 0 ) t + ε 2 μ N + ε 2 V 2 ( t ) V 2 ( 0 ) t + δ μ N + δ R ( t ) R ( 0 ) t = b ( μ N + q 0 ) S t k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t + ( q 0 ε 1 ( μ N + ε 2 ) + ε 2 q 1 ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ) S t .
Then,
S t = 1 ( μ N + q 0 ) q 0 ε 1 ( μ N + ε 2 ) + ε 2 q 1 ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ) φ ( t ) = b μ N ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) + ( μ N q 0 + q 0 q 1 + ε 2 q 0 ) 1 μ N ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) + ( μ N q 0 + q 0 q 1 + ε 2 q 0 ) k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t φ ( t ) = S 0 S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t φ ( t ) ,  
where
φ ( t ) = S 0 b [ ( S ( t ) S ( 0 ) t + k = 1 K I k ( t ) I k ( 0 ) t ) + ε 1 ( μ N + ε 2 ) + ε 2 q 1 ( μ N + ε 1 + q 1 ) ( μ N + ε 2 ) V 1 ( t ) V 1 ( 0 ) t + ε 2 μ N + ε 2 V 2 ( t ) V 2 ( 0 ) t + δ μ N + δ R ( t ) R ( 0 ) t ] .
Note that lim t N ( t ) t lim t 1 t b μ N = 0 , we have
lim t φ ( t ) = 0 ,   a . s .
Applying It o ^ ’s formula to the equations of I k in model (12) and using (16) again, we obtain
ln I k ( t ) ln I k ( 0 ) t = β 0 k S t ( μ N + μ k + γ k ) σ k 2 2 S 2 t + 1 t σ k 0 t S d B k ( s ) ,   a . s .
By the Cauchy–Schwarz inequality, S 2 t S t 2   a . s . Then,
ln I k ( t ) ln I k ( 0 ) t β 0 k S t ( μ N + μ k + γ k ) σ k 2 2 S t 2 + 1 t σ k 0 t S d B k ( s ) ,   a . s .
Now, using (A12) and (A16),
ln I k ( t ) ln I k ( 0 ) t β 0 k [ S 0 S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t φ ( t ) ] ( μ N + μ k + γ k ) σ k 2 2 [ ( S 0 ) 2 + ( S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ) 2 + φ 2 ( t ) 2 S 0 b ( S 0 φ ( t ) ) k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t 2 S 0 φ ( t ) ] + 1 t σ k 0 t S d B k ( s )
We define the last term of Equation (A17) as
M k ( t ) = σ k 0 t S d B k ( s ) .
Then, M k , M k t = 0 t σ k 2 S 2 ( r ) d r σ k 2 ( b μ N ) 2 t , and we have
limsup t M k , M k t t σ k 2 ( b μ N ) 2 < .
By Lemma 1,
lim t M k ( t ) t = 0 ,   a . s .
Plugging M k ( t ) into Equation (A17),
ln I k ( t ) t ln I k ( 0 ) t + β 0 k [ S 0 S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t φ ( t ) ] ( μ N + μ k + γ k ) σ k 2 2 [ S 02 + ( S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ) 2 + φ 2 ( t ) 2 S 0 b ( S 0 φ ( t ) ) k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t 2 S 0 φ ( t ) ] + M k ( t ) t ln I k ( 0 ) t + β 0 k [ S 0 φ ( t ) ] ( μ N + μ k + γ k ) σ k 2 2 [ S 02 + ( S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ) 2 + φ 2 ( t ) + 2 φ ( t ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t 2 S 0 φ ( t ) ] ( β 0 k S 0 σ k 2 ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t + M k ( t ) t .
By the assumption β 0 k > S 0 σ k 2 , we have
ln I k ( t ) t ln I k ( 0 ) t + β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 S 02 β 0 k φ ( t ) σ k 2 2 [ 2 φ ( t ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t 2 S 0 φ ( t ) ] .
As R 0 k S 1 = R 0 k σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 1 < 0 , we have
limsup t ln I k ( t ) t β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( S 0 ) 2 = ( μ N + μ k + γ k ) [ R 0 k 1 σ k 2 2 ( μ N + μ k + γ k ) ( S 0 ) 2 ] < 0 ,   a . s .
(b)
If β 0 k 2 2 ( μ N + μ k + γ k ) < σ k 2 , by Equation (A16),
ln I k ( t ) ln I k ( 0 ) t β 0 k S t ( μ N + μ k + γ k ) σ k 2 2 S t 2 + 1 t σ k 0 t S d B k ( s ) = ( μ N + μ k + γ k ) σ k 2 2 ( S t β 0 k σ k 2 ) 2 + β 0 k 2 2 σ k 2 + 1 t σ k 0 t S d B k ( s ) ( μ N + μ k + γ k ) + β 0 k 2 2 σ k 2 + 1 t σ k 0 t S d B k ( s ) .
It implies
limsup t ln I k ( t ) t ( μ N + μ k + γ k ) + β 0 k 2 2 σ k 2 < 0 ,   a . s .
The proof is complete. □

Appendix C

Proof of Theorem 4.
(a)
Using (A12) and (A16),
ln I k ( t ) ln I k ( 0 ) t = β 0 k S t ( μ N + μ k + γ k ) σ k 2 2 S 2 t + 1 t σ k 0 t S d B k ( s ) β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( b μ N ) 2 β 0 k S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t + β 0 k φ ( t ) + 1 t σ k 0 t S d B k ( s ) ,
We have
β 0 k S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ln I k ( t ) ln I k ( 0 ) t + β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( b μ N ) 2 + β 0 k φ ( t ) + 1 t σ k 0 t S d B k ( s ) ,
k = 1 , 2 , , N . An upper bound of the left side of inequality (A27) is
β 0 k S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) b μ N ,
implying limsup t ( ln I k ( t ) ln I k ( 0 ) t ) has an upper bound. As I k ( t ) ( 0 , b μ N ) , we have
liminf t ( ln I k ( t ) ln I k ( 0 ) t ) 0 ,   a . s .
Then,
liminf t k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ( μ N + μ k + γ k ) ( R 0 k * 1 ) β 0 k S 0 b ,
k = 1 , 2 , , N . We set k ^ * = arg max k : R 0 k * > 1 { ( μ N + μ k + γ k ) ( R 0 k * 1 ) β 0 k } . Then
liminf t k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ( μ N + μ k ^ * + γ k ^ * ) ( R 0 k ^ * * 1 ) β 0 k ^ * S 0 b ,   a . s .
Setting k ˜ * = arg max k : R 0 k * > 1 { ( μ N + μ k + γ k ) δ γ k μ N + δ } , Equation (22) is held.
(b)
Using (A12) and (A17),
ln I k ( t ) ln I k ( 0 ) t β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( S 0 ) 2 + ϕ ( t ) ( β 0 k S 0 σ k 2 ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t + M k ( t ) t ,
where
ϕ ( t ) = β 0 k φ ( t ) σ k 2 2 [ φ 2 ( t ) + 2 φ ( t ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t 2 S 0 φ ( t ) ] .
By Equation (A14), lim t φ ( t ) = 0 implies lim t ϕ ( t ) = 0 ,   a . s . Then, we have
( β 0 k S 0 σ k 2 ) S 0 b k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ln I k ( t ) ln I k ( 0 ) t + β 0 k S 0 ( μ N + μ k + γ k ) σ k 2 2 ( S 0 ) 2 + ϕ ( t ) + M k ( t ) t .  
By assumption, there exists some k 0 { 1 , 2 , , K } such that liminf t I k 0 t > 0 , β 0 k 0 > S 0 σ k 0 2 , and R 0 k 0 S R 0 k 0 * > 1 . With I k 0 b μ N , we obtain lim t ln I k 0 ( t ) ln I k 0 ( 0 ) t = 0 ,   a . s . Thus,
( β 0 k 0 S 0 σ k 0 2 ) S 0 b limsup t k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t β 0 k 0 S 0 ( μ N + μ k 0 + γ k 0 ) σ k 2 2 ( S 0 ) 2 .
We have
limsup t k = 1 K ( ( μ N + μ k + γ k ) δ γ k μ N + δ ) I k t ( μ N + μ k ^ 0 + γ k ^ 0 ) ( R 0 k ^ 0 S 1 ) ( β 0 k ^ 0 S 0 σ k ^ 0 2 ) S 0 b ,
where k ^ 0 : = argmax k 0 : R 0 k 0 S > 1 { ( μ N + μ k 0 + γ k 0 ) ( R 0 k 0 S 1 ) ( β 0 k 0 S 0 σ k 0 2 )   } .
By the definition of k ˜ * , we have
limsup t ( ( μ N + μ k ˜ * + γ k ˜ * ) δ γ k ˜ * μ N + δ ) k = 1 K I k t ( μ N + μ k ^ 0 + γ k ^ 0 ) ( R 0 k ^ 0 S 1 ) ( β 0 k ^ 0 S 0 σ k ^ 0 2 ) S 0 b
and then Equation (23) is held.
The proof is complete. □

Appendix D

Proof of Theorem 5.
Using Equation (A15),
β 0 k ¯ * ln I k ( t ) ln I k ( 0 ) t β 0 k ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t                 = ( β 0 k ¯ * β 0 k β 0 k β 0 k ¯ * ) S t β 0 k ¯ * ( μ N + μ k + γ k ) + β 0 k ( μ N + μ k ¯ * + γ k ¯ * )                                               ( β 0 k ¯ * σ k 2 2 β 0 k σ k ¯ * 2 2 ) S 2 t + 1 t β 0 k ¯ * σ k 0 t S d B k ( s ) 1 t β 0 k σ k ¯ * 0 t S d B k ¯ * ( s )                   = ( μ N + μ k + γ k ) ( μ N + μ k ¯ * + γ k ¯ * ) ×                         [ β 0 k μ N + μ k + γ k β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * ( β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 2 ( μ N + μ k ¯ * + γ k ¯ * ) ) S 2 t ] + 1 t β 0 k ¯ * σ k 0 t S d B k ( s ) 1 t β 0 k σ k ¯ * 0 t S d B k ¯ * ( s )
By assumptions,
lim t ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t = lim t ( β 0 k ¯ * S t ( μ N + μ k ¯ * + γ k ¯ * ) σ k ¯ * 2 2 S 2 t + 1 t σ k ¯ * 0 t S d B k ¯ * )                                                             = ( μ N + μ k ¯ * + γ k ¯ * ) lim t ( β 0 k ¯ * S t μ N + μ k ¯ * + γ k ¯ * σ k ¯ * 2 2 ( μ N + μ k ¯ * + γ k ¯ * ) S 2 t 1 ) = 0 ,   a . s .
There exists a large enough T 0 > 0 such that for all t T 0 , and all ω Ω 0 Ω with P ( Ω 0 ) = 1 , we have
S 2 t ( ω ) = 2 β 0 k ¯ * S t ( ω ) σ k ¯ * 2 2 ( μ N + μ k ¯ * + γ k ¯ * ) σ k ¯ * 2 + ϵ t ( ω ) ,
where | ϵ t ( ω ) | 0 and | ϵ t ( ω ) | 0 when t .
Substituting (A39) to the first term on the right-hand side of (A37), we have
β 0 k ¯ * ln I k ( t ) ln I k ( 0 ) t β 0 k ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t = β 0 k ¯ * ( μ N + μ k + γ k ) ( μ N + μ k ¯ * + γ k ¯ * ) σ k ¯ * 2 [ σ k ¯ * 2 μ N + μ k ¯ * + γ k ¯ * ( β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) ) S t ( ω ) + σ k 2 ( μ N + μ k + γ k ) ] + ϵ t * ( ω ) ,
where
ϵ t * ( ω ) = ϵ t ( ω ) ( β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 2 ( μ N + μ k ¯ * + γ k ¯ * ) )                                                                                                     + 1 t β 0 k ¯ * σ k 0 t S d B k ( s , ω ) 1 t β 0 k σ k ¯ * 0 t S d B k ¯ * ( s , ω ) .
Now, if β 0 k ¯ * σ k ¯ * 2 β 0 k σ k 2 , then β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) > 0 . As S t 0   a . s . and σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) > σ k 2 ( μ N + μ k + γ k ) , the right hand side of Equation (A40) is negative when ϵ t * ( ω ) is dropped. Then,
limsup t β 0 k ¯ * ln I k ( t ) ln I k ( 0 ) t β 0 k ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t < 0   a . s .
For lim t ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t = 0 ,   a . s . , we have β 0 k ¯ * limsup t ln I k ( t ) ln I k ( 0 ) t < 0 ,   a . s .
If β 0 k ¯ * σ k ¯ * 2 < β 0 k σ k 2 , then β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) < 0 . By Equation (A12), there exists a T 1 > 0 large enough such that S t S 0 ,   a . s .   t T 1 . Thus, we have Ω 1 Ω , P ( Ω 1 ) = 1 , such that for all ω Ω 0 Ω 1 and t > max { T 0 , T 1 } ,
β 0 k ¯ * ln I k ( t ) ln I k ( 0 ) t β 0 k ln I k ¯ * ( t ) ln I k ¯ * ( 0 ) t β 0 k ¯ * ( μ N + μ k + γ k ) ( μ N + μ k ¯ * + γ k ¯ * ) σ k ¯ * 2 [ σ k ¯ * 2 μ N + μ k ¯ * + γ k ¯ * ( β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 ( μ N + μ k + γ k ) β 0 k μ N + μ k + γ k σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) ) S 0 + σ k 2 ( μ N + μ k + γ k ) ] + ϵ t * ( ω ) .
The part in the rectangle bracket on the right-hand side of Equation (A43) is
σ k ¯ * 2 μ N + μ k ¯ * + γ k ¯ * + ( β 0 k μ N + μ k + γ k σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) β 0 k ¯ * μ N + μ k ¯ * + γ k ¯ * σ k 2 ( μ N + μ k + γ k ) ) S 0 σ k 2 ( μ N + μ k + γ k ) = ( β 0 k S 0 μ N + μ k + γ k 1 ) σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) ( β 0 k ¯ * S 0 μ N + μ k ¯ * + γ k ¯ * 1 ) σ k 2 ( μ N + μ k + γ k ) = ( R 0 k S 1 ) σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) ( R 0 k ¯ * S 1 ) σ k 2 ( μ N + μ k + γ k ) .
Then, if σ k 2 ( μ N + μ k + γ k ) > ( R 0 k S 1 R 0 k ¯ * S 1 ) σ k ¯ * 2 ( μ N + μ k ¯ * + γ k ¯ * ) , the right hand side of Equation (A43) is negative when ϵ t * ( ω ) is dropped.
Again, it implies limsup t ln I k ( t ) ln I k ( 0 ) t < 0 ,   a . s .
The proof is complete. □

References

  1. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A Contain. Pap. A Math. Phys. Character 1927, 115, 700–721. [Google Scholar]
  2. Anderson, R.M.; May, R.M. Population biology of infectious diseases: Part I. Nature 1979, 280, 361–367. [Google Scholar] [CrossRef] [PubMed]
  3. May, R.M.; Anderson, R.M. Population biology of infectious diseases: Part II. Nature 1979, 280, 455–461. [Google Scholar] [CrossRef] [PubMed]
  4. Capasso, V.; Serio, G. A generalization of the Kermack-McKendrick deterministic epidemic model. Math. Biosci. 1978, 42, 43–61. [Google Scholar] [CrossRef]
  5. Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: Boston, MA, USA, 2015. [Google Scholar]
  6. Gray, A.; Greenhalgh, D.; Hu, L.; Mao, X.; Pan, J. A stochastic differential equation SIS epidemic model. SIAM J. Appl. Math. 2011, 71, 876–902. [Google Scholar] [CrossRef] [Green Version]
  7. Méndez, V.; Campos, D.; Horsthemke, W. Stochastic fluctuations of the transmission rate in the susceptible-infected-susceptible epidemic model. Phys. Rev. E 2012, 86, 011919. [Google Scholar] [CrossRef] [Green Version]
  8. Cai, Y.; Kang, Y.; Banerjee, M.; Wang, A. A stochastic SIRS epidemic model with infectious force under intervention strategies. J. Differ. Equ. 2015, 259, 7463–7502. [Google Scholar] [CrossRef]
  9. Chen, C.; Kang, Y. Dynamics of a stochastic multi-strain SIS epidemic model driven by Levy noise. Commun. Nonlinear Sci. Numer. Simul. 2017, 43, 379–395. [Google Scholar] [CrossRef]
  10. Hethcote, H.; Ma, Z.; Liao, S. Effects of quarantine in six endemic models for infectious diseases. Math. Biosci. 2002, 180, 141–160. [Google Scholar] [CrossRef]
  11. Huang, S.; Chen, F.; Chen, L. Global dynamics of a network-based SIQRS epidemic model with demographics and vaccination. Commun. Nonlinear Sci. Numer. Simul. 2017, 43, 296–310. [Google Scholar] [CrossRef]
  12. Li, K.; Zhu, G.; Ma, Z.; Chen, L. Dynamic stability of an SIQS epidemic network and its optimal control. Commun. Nonlinear Sci. Numer. Simul. 2019, 66, 84–95. [Google Scholar] [CrossRef]
  13. Zhao, R.; Liu, Q.; Sun, M. Dynamical behavior of a stochastic SIQS epidemic model on scale-free networks. J. Appl. Math. Comput. 2022, 68, 813–838. [Google Scholar] [CrossRef]
  14. Zhang, X.B.; Zhang, X.H. The threshold of a deterministic and stochastic SIQS epidemic model with varying total population size. Appl. Math. Model. 2021, 91, 749–767. [Google Scholar] [CrossRef]
  15. Li, J.; Ma, Z. Qualitative analyses of SIS epidemic model with vaccination and varying total population size. Math. Comput. Model. 2002, 35, 1235–1243. [Google Scholar] [CrossRef]
  16. Li, J.; Ma, Z. Global analysis of SIS epidemic models with variable total population size. Math. Comput. Model. 2004, 39, 1231–1242. [Google Scholar]
  17. Zhao, Y.; Zhang, Q.; Jiang, D. The asymptotic behavior of a stochastic SIS epidemic model with vaccination. Adv. Differ. Equ. 2015, 2015, 328. [Google Scholar] [CrossRef] [Green Version]
  18. Zaman, G.; Kang, Y.H.; Jung, I.H. Stability analysis and optimal vaccination of an SIR epidemic model. Biosystems 2008, 93, 240–249. [Google Scholar] [CrossRef]
  19. Fudolig, M.; Howard, R. The local stability of a modified multi-strain sir model for emerging viral strains. PLoS ONE 2020, 15, e0243408. [Google Scholar] [CrossRef]
  20. World Health Organization. Coronavirus Disease (COVID-2019) Situation Reports. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports (accessed on 22 April 2022).
  21. Dong, Y.; Dai, T.; Wei, Y.; Zhang, L.; Zheng, M.; Zhou, F. A systematic review of SARS-CoV-2 vaccine candidates. Signal Transduct. Target. Ther. 2020, 5, 237. [Google Scholar] [CrossRef]
  22. Rella, S.A.; Kulikova, Y.A.; Dermitzakis, E.T.; Kondrashov, F.A. Rates of SARS-CoV-2 transmission and vaccination impact the fate of vaccine-resistant strains. Sci. Rep. 2021, 11, 15729. [Google Scholar] [CrossRef]
  23. Zheng, Q.; Wang, X.; Bao, C.; Ma, Z.; Pan, Q. Mathematical modelling and projecting the second wave of COVID-19 pandemic in Europe. J. Epidemiol. Community Health 2021, 75, 601–603. [Google Scholar] [CrossRef] [PubMed]
  24. Zheng, Q.; Wang, X.; Bao, C.; Ji, Y.; Liu, H.; Meng, Q.; Pan, Q. A multi-regional hierarchical-tier mathematical model of the spread and control of COVID-19 epidemics from epicentre to adjacent regions. Transbound. Emerg. Dis. 2022, 69, 549–558. [Google Scholar] [CrossRef] [PubMed]
  25. Zheng, Q.; Wang, X.; Pan, Q.; Wang, L. Optimal strategy for a dose-escalation vaccination against COVID-19 in refugee camps. AIMS Math. 2022, 7, 9288–9310. [Google Scholar] [CrossRef]
  26. Blower, S.M.; Dowlatabadi, H. Sensitivity and uncertainty analysis of complex models of disease transmission: An HIV model, as an example. Int. Stat. Rev. 1994, 62, 229–243. [Google Scholar] [CrossRef]
  27. Zhang, K.; Wang, X.; Liu, H.; Ji, Y.; Pan, Q.; Wei, Y.; Ma, M. Mathematical analysis of a human papillomavirus transmission model with vaccination and screening. AIMS Math. Biosci. Eng. 2020, 17, 5449–5476. [Google Scholar] [CrossRef]
  28. Husnulkhotimah, H.; Rusin, R.; Aldila, D. Using the SVIRS model to understand the prevention strategy for influenza with vaccination. AIP Conf. Proc. 2021, 2374, 030004. [Google Scholar]
  29. Mao, X. Stochastic Differential Equations and Applications; Horwood Publishing: Chichester, UK, 1997; pp. 12–13. [Google Scholar]
Figure 1. Compartment diagram for the proposed multi-strain SIR model.
Figure 1. Compartment diagram for the proposed multi-strain SIR model.
Mathematics 10 01804 g001
Figure 2. Significance test of 3-train model parameters and PRCC results for R 0 .
Figure 2. Significance test of 3-train model parameters and PRCC results for R 0 .
Mathematics 10 01804 g002
Figure 3. The simulated numbers of infected individuals under condition (a) of Theorem 2.
Figure 3. The simulated numbers of infected individuals under condition (a) of Theorem 2.
Mathematics 10 01804 g003
Figure 4. The simulated numbers of infected individuals under condition (b) of Theorem 2.
Figure 4. The simulated numbers of infected individuals under condition (b) of Theorem 2.
Mathematics 10 01804 g004
Figure 5. The simulated numbers of infected individuals under the conditions of the disease persistence.
Figure 5. The simulated numbers of infected individuals under the conditions of the disease persistence.
Mathematics 10 01804 g005
Figure 6. The simulated numbers of total infected individuals under the different ( q 0 , q 1 ) combinations.
Figure 6. The simulated numbers of total infected individuals under the different ( q 0 , q 1 ) combinations.
Mathematics 10 01804 g006
Figure 7. The simulated numbers of infected individuals under the conditions of Table 1 in [22].
Figure 7. The simulated numbers of infected individuals under the conditions of Table 1 in [22].
Mathematics 10 01804 g007
Table 1. Description of variables and parameters in model (1).
Table 1. Description of variables and parameters in model (1).
Variable/ParameterDescription
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   i th   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
Table 2. Partial rank correlation coefficients.
Table 2. Partial rank correlation coefficients.
ParameterPRCCParameterPRCC
b 0.8258 ** γ 2 −0.1198 **
β 1 0.3185 ** γ 3 −0.1515 **
β 2 0.2817 ** μ 1 −0.2133 **
β 3 0.2523 ** μ 2 −0.1224 **
δ −0.0699 * μ 3 −0.0917 **
ε 1 0.0530 μ N −0.8421 **
ε 2 0.1044 ** q 0 −0.3632 **
γ 1 −0.1843 ** q 1 −0.0481
The results are significant at the 0.05 level (*) or the 0.01 level (**).
Table 3. The different ( q 0 ,   q 1 ) combinations and related R 0 S ’s and ξ ¯ ’s.
Table 3. The different ( q 0 ,   q 1 ) combinations and related R 0 S ’s and ξ ¯ ’s.
( q 0 , q 1 ) R 0 S ξ ¯
(0.1, 0.1)8.66625.7327
(0.1, 0.5)8.13005.6695
(0.1, 0.9)7.94105.6458
(0.5, 0.1)4.58334.9769
(0.5, 0.5)3.90224.7210
(0.5, 0.9)3.69124.6240
(0.9, 0.1)3.11514.2961
(0.9, 0.5)2.56693.8537
(0.9, 0.9)2.40413.6848
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chang, Y.-C.; Liu, C.-T. A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate. Mathematics 2022, 10, 1804. https://doi.org/10.3390/math10111804

AMA Style

Chang Y-C, Liu C-T. A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate. Mathematics. 2022; 10(11):1804. https://doi.org/10.3390/math10111804

Chicago/Turabian Style

Chang, Yen-Chang, and Ching-Ti Liu. 2022. "A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate" Mathematics 10, no. 11: 1804. https://doi.org/10.3390/math10111804

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop