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

Chang, Yen-Chang; Liu, Ching-Ti

doi:10.3390/math10111804

Open AccessArticle

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

by

Yen-Chang Chang

^1,*,† and

Ching-Ti Liu

^2,*,†

¹

Center for General Education, National Tsing Hua University, Hsinchu City 300, Taiwan

²

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

(This article belongs to the Special Issue Recent Advances of Computational Statistics in Industry and Business II)

Download

Browse Figures

Versions Notes

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.

Keywords:

SIR model; vaccination; basic reproduction numbers; extinction; persistence

MSC:

60H35; 37A99

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

{\begin{cases} d S (t) = [b - S (t) \sum_{k = 1}^{K} β_{k} I_{k} (t) - (μ_{N} + q_{0}) S (t) + \sum_{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, \dots, K, \\ d R (t) = [\sum_{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 . \end{cases}

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

Suppose that the number of total populations is

N (t) = S (t) + \sum_{k = 1}^{K} I_{k} (t) + \sum_{i = 1}^{2} V_{i} (t) + R (t)

. It implies

\begin{array}{l} d [S (t) + \sum_{k = 1}^{K} I_{k} (t) + R (t) + \sum_{i = 1}^{2} V_{i} (t)] \\ = [b - μ_{N} N (t) - \sum_{k = 1}^{K} μ_{k} I_{k} (t)] d t \leq [b - μ_{N} N (t)] d t \end{array}

(2)

Then, we have

N (t) \leq N_{1}

,

\forall t \geq 0

, where

N_{1} = {\begin{matrix} \frac{b}{μ_{N}}, if N (0) \leq \frac{b}{μ_{N}}; \\ N (0), if N (0) > \frac{b}{μ_{N}} . \end{matrix}

(3)

To simplify the following discussion, we assume

N (0) \leq \frac{b}{μ_{N}}

. Next, we will discuss the stability of model (1) on an invariant set

Λ

:

Λ = {θ = (θ_{1}, \dots, θ_{K + 4}) \in ℝ_{+}^{K + 4} | \sum_{k = 1}^{K + 4} θ_{k} \leq \frac{b}{μ_{N}}},

(4)

where

θ_{1} = S

,

(θ_{2}, \dots, θ_{K + 1}) = I = (I_{1}, \dots, 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}, \dots, I_{k}^{0}, R^{0}, V_{1}^{0}, V_{2}^{0})

, where

{\begin{array}{l} S^{0} = \frac{b}{μ_{N}} \frac{(μ_{N} + ε_{2}) (μ_{N} + ε_{1} + q_{1})}{(μ_{N} + ε_{2}) (μ_{N} + ε_{1} + q_{1}) + q_{0} (μ_{N} + ε_{2}) + q_{0} q_{1}}, \\ V_{1}^{0} = \frac{b}{μ_{N}} \frac{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} = \frac{b}{μ_{N}} \frac{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} = \dots = I_{K}^{0} = R^{0} = 0 . \end{array}

(5)

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^{*}, {\tilde{I}}_{k}^{*}, R^{*}, V_{1}^{*}, V_{2}^{*})

, and

{\begin{array}{l} S^{*} = \frac{μ_{N} + μ_{k} + γ_{k}}{β_{k}}, \\ V_{1}^{*} = \frac{q_{0} (μ_{N} + μ_{k} + γ_{k})}{β_{0 k} (μ_{N} + ε_{1} + q_{1})}, \\ V_{2}^{*} = \frac{q_{0} q_{1} (μ_{N} + μ_{k} + γ_{k})}{β_{k} (μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}, \\ R^{*} = \frac{γ_{k} I_{k}^{*}}{μ_{N} + δ}, \\ I_{k}^{*} = \frac{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} + δ))} \end{array}

(6)

k = 1, 2, \dots, K

. Note that

{\tilde{I}}_{k}^{*} = (I_{1}^{*}, \dots, I_{K}^{*})

means

I_{k}^{*} > 0

and

I_{l}^{*} = 0, \forall l \neq k

.

Using the next generation method, we can obtain the basic reproduction number.

Set

F = diag (β_{1} S, \dots β_{K} S)

(7)

and

V = diag (μ_{N} + μ_{1} + γ_{1}, \dots μ_{N} + μ_{K} + γ_{K})

(8)

The next generation matrix is

F V^{- 1} = diag (\frac{β_{k} S}{μ_{N} + μ_{k} + γ_{k}}, k = 1, 2, \dots, K)

(9)

Substituting

S

with

S^{0}

, we have the basic reproduction number

R_{0 k} = \frac{β_{k} S^{0}}{μ_{N} + μ_{k} + γ_{k}}

for strain

k

in the deterministic model, that is,

R_{0 k} : = \frac{b}{μ_{N}} \frac{β_{k} (μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}{(μ_{N} + μ_{k} + γ_{k}) [(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2}) + q_{0} (μ_{N} + ε_{2}) + q_{0} q_{1}]},

(10)

k = 1, 2, \dots, K

. This yields the basic reproduction number of model (1) with

R_{0} = \max {R_{0 k}, k = 1, 2, \dots, K}

(11)

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 \geq 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, \dots, K

. Thus, the proposed stochastic model becomes

{\begin{cases} d S (t) = [b - S \sum_{k = 1}^{K} β_{0 k} I_{k} - (μ_{N} + q_{0}) S + \sum_{i = 1}^{2} ε_{i} V_{i} + δ R] d t - S \sum_{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, \dots, K, \\ d R (t) = [\sum_{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 . \end{cases}

(12)

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) \in Λ

, model (12) has a unique positive solution

θ (t) \in Λ

for all

t \geq 0

with probability 1. That is,

P (θ (t) \in Λ, \forall t \geq 0) = 1 .

(13)

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} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2},

(14)

k = 1, 2, \dots, K

. The basic reproductive number for model (12) is

R_{0}^{S} : = \max {R_{0 k}^{S}, k = 1, 2, \dots, K} .

(15)

To facilitate subsequent discussions, we define

x_{t} = {〈 x 〉}_{t} = \frac{1}{t} \int_{0}^{t} x (s) d s and {〈 x, x 〉}_{t} = \frac{1}{t} \int_{0}^{t} x^{2} (s) d s .

(16)

The following lemma is useful for our discussion.

Lemma 1

(Strong Law of Large Numbers [29]). Let

M = {M_{t}}_{t \geq 0}

be a real-valued continuous local martingale vanishing at

t = 0

. Then

\lim_{t \to \infty} {〈 M, M 〉}_{t} = \infty, a . s . \Rightarrow \lim_{t \to \infty} \frac{M_{t}}{{〈 M, M 〉}_{t}} = 0, a . s .

(17)

and also

\underset{t \to \infty}{limsup} \frac{{〈 M, M 〉}_{t}}{t} < \infty, a . s . \Rightarrow \lim_{t \to \infty} \frac{M_{t}}{t} = 0, a . s .

(18)

Theorem 2.

In model (12), consider the following two conditions:

(a): $R_{0 k}^{S} = R_{0 k} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2} < 1$ and $β_{0 k} > S^{0} σ_{k}^{2}$ , or
(b): $\frac{β_{0 k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} < σ_{k}^{2}$ .

If one of the above conditions holds, then for an arbitrary initial vector

θ (0) \in Λ

, the solution

θ (t) \in Λ

obeys

\underset{t \to \infty}{limsup} \frac{\ln I_{k} (t)}{t} < 0, a . s .

(19)

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, \dots, 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} = \frac{β_{0 k} S^{0}}{(μ_{N} + μ_{k} + γ_{k})} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2} .

(20)

Note that

f (\frac{β_{0 k}}{σ_{k}^{2}}) = \frac{β_{0 k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k}) σ_{k}^{2}}

is the maximum of

f

. Then,

\frac{β_{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} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2} < 1$ and $β_{0 k} > S^{0} σ_{k}^{2}$ , or
(b): $\frac{β_{0 k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} < σ_{k}^{2}$ ,

If one of the above conditions holds,then for arbitrary initial vector

θ (0) \in Λ

, the solution

θ (t) \in Λ

obeys

\underset{t \to \infty}{limsup} \frac{\ln I_{k} (t)}{t} < 0, a . s .

That is,

I_{k}

tends to zero exponentially with probability 1 for

k = 1, 2, \dots, K

.

3.3. Persistence

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

\underset{t \to \infty}{liminf} \sum_{k = 1}^{K} I_{k} > 0, a . s .

holds. We define

R_{0 k}^{*} : = R_{0 k} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(\frac{b}{μ_{N}})}^{2},

(21)

k = 1, 2, \dots, K

. As

R_{0 k}^{*} \leq R_{0 k}^{S}

,

R_{0 k}^{*} > 1

implies

R_{0 k}^{S} > 1

.

Theorem 4.

For arbitrary initial vector

θ (0) \in Λ

,

(a): If $\max {R_{0 k}^{*}} > 1$ , we have

$\underset{t \to \infty}{liminf} \sum_{k = 1}^{K} I_{k}_{t} \geq \underline{ξ} : = \frac{(μ_{N} + μ_{{\hat{k}}^{*}} + γ_{{\hat{k}}^{*}}) (R_{0 {\hat{k}}^{*}}^{*} - 1)}{β_{0 {\hat{k}}^{*}} \frac{S^{0}}{b} ((μ_{N} + μ_{{\tilde{k}}^{*}} + γ_{{\tilde{k}}^{*}}) - \frac{δ γ_{{\tilde{k}}^{*}}}{μ_{N} + δ})} a . s .,$

(22)

where ${\tilde{k}}^{*} = \underset{k : R_{0 k}^{*} > 1}{\arg \max} {(μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}}$ and ${\hat{k}}^{*} = \underset{k : R_{0 k}^{*} > 1}{\arg \max} {\frac{(μ_{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

$\underset{t \to \infty}{limsup} \sum_{k = 1}^{K} 〈 I_{k} 〉_{t} \leq \bar{ξ} : = \frac{(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) (R_{0 {\hat{k}}_{0}}^{S} - 1)}{(β_{0 {\hat{k}}_{0}} - S^{0} σ_{{\hat{k}}_{0}}^{2}) \frac{S^{0}}{b} ((μ_{N} + μ_{{\tilde{k}}_{*}} + γ_{{\tilde{k}}_{*}}) - \frac{δ γ_{{\tilde{k}}_{*}}}{μ_{N} + δ})} a . s .,$

(23)

where ${\tilde{k}}_{*} = \underset{k : R_{0 k}^{*} > 1}{\arg \min} {(μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}}$ and ${\hat{k}}_{0} = \underset{k_{0} : R_{0 k_{0}}^{S} > 1}{a r g m a x} {\frac{(μ_{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

\sum_{k = 1}^{K} 〈 I_{k} 〉_{t}

should oscillate in the interval

[\underline{ξ}, \bar{ξ}]

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

{\bar{k}}^{*} = \arg \max_{k : R_{0 k}^{S} > 1} R_{0 k}^{S} .

(24)

Then, it is obvious that

R_{0 {\bar{k}}^{*}}^{S} = R_{0}^{S}

. We prove that, under specified conditions, strain

{\bar{k}}^{*}

is the uniquely persistent strain in Theorem 4.

Theorem 5.

Suppose the conditions of Theorem 4 all hold. Assume that

\lim_{t \to \infty} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} = 0, a . s .

(25)

and there are some

k \neq {\bar{k}}^{*}

with

R_{0 {\bar{k}}^{*}}^{S} > R_{0 k}^{S} > 1

. If

(a): $\frac{β_{0 {\bar{k}}^{*}}}{σ_{{\bar{k}}^{*}}^{2}} \geq \frac{β_{0 k}}{σ_{k}^{2}}$ and $\frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} > \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})}$ , or
(b): $\frac{β_{0 {\bar{k}}^{*}}}{σ_{{\bar{k}}^{*}}^{2}} < \frac{β_{0 k}}{σ_{k}^{2}}$ and $\frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} > (\frac{R_{0 k}^{S} - 1}{R_{0 {\bar{k}}^{*}}^{S} - 1}) \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}$ ,

we have

\underset{t \to \infty}{limsup} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} < 0, a . s .

, i.e.,

\lim_{t \to \infty} 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 = \frac{b}{μ_{N}}, B = μ_{N}^{2} + μ_{N} (ε_{2} + ε_{1}) + ε_{1} ε_{2}

, and

C = ε_{2} + μ_{N}

, then

S^{0} = A \frac{B + q_{1} (ε_{2} + μ_{N})}{B + q_{1} (ε_{2} + μ_{N}) + q_{0} (ε_{2} + μ_{N}) + q_{0} q_{1}} = A \frac{B^{'} + q_{1}}{B^{'} + q_{1} + q_{0} + q_{1} q_{0} / C},

(26)

where

B^{'} = \frac{B}{C} = \frac{μ_{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 \frac{B^{'} + q_{1}}{B^{'} + q_{1} + q_{1}^{2} + q_{1}^{3} / C}

, and

\frac{\partial S^{0}}{\partial q_{1}} = A (\frac{1}{B^{'} + q_{1} + q_{1}^{2} + \frac{a q_{1}^{3}}{C}} - \frac{(B^{'} + q_{1}) (1 + 2 q_{1} + 3 \frac{q_{1}^{2}}{C})}{{(B^{'} + q_{1} + q_{1}^{2} + \frac{q_{1}^{3}}{C})}^{2}}) = - A (\frac{2 B^{'} q_{1} + (\frac{3 B^{'}}{C} + 1) q_{1}^{2} + \frac{2 q_{1}^{3}}{C}}{{(B^{'} + q_{1} + q_{1}^{2} + \frac{q_{1}^{3}}{C})}^{2}}) < 0 .

(27)

As discussed in Remark 1,

R_{0 k}^{S}

has a global maximum at

S^{0} = \frac{β_{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

\frac{\partial R_{0 k}^{S}}{\partial q_{1}} = \frac{\partial R_{0 k}^{S}}{\partial S^{0}} \frac{\partial S^{0}}{\partial q_{1}} = \frac{- A}{μ_{N} + μ_{k} + γ_{k}} (β_{k} - σ_{k}^{2} S^{0}) (\frac{2 B^{'} q_{1} + (\frac{3 B^{'}}{C} + 1) q_{1}^{2} + \frac{2 q_{1}^{3}}{C}}{{(B^{'} + q_{1} + q_{1}^{2} + q_{1}^{3} / C)}^{2}}) < 0 .

(28)

R_{0 k}^{S}

is a decreasing function of

q_{0}

. For

q_{1} \in (0, 1)

, we have

S^{0} \in (\frac{b}{μ_{N}} \frac{μ_{N} + ε_{1} + 1}{μ_{N} + ε_{1} + 2 + 1 / (ε_{2} + μ_{N})}, \frac{b}{μ_{N}})

and

R_{0 k}^{S} \in \frac{A}{μ_{N} + μ_{k} + γ_{k}} \times (β_{0 k} (\frac{B^{'} + 1}{B^{'} + 2 + \frac{1}{C}}) - \frac{σ_{k}^{2}}{2} A {(\frac{B^{'} + 1}{B^{'} + 2 + \frac{1}{C}})}^{2}, β_{0 k} - \frac{σ_{k}^{2}}{2} A), k = 1, 2, \dots, K .

(29)

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}),

(30)

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}) \leq C_{0} and 0 \leq q_{1} \leq 1

(31)

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 (

\bar{ξ}

) as an estimate of the number of infected individuals. Let

D = \frac{1}{b} ((μ_{N} + μ_{{\tilde{k}}_{*}} + γ_{{\tilde{k}}_{*}}) - \frac{δ γ_{{\tilde{k}}_{*}}}{μ_{N} + δ}) .

(32)

Then,

\bar{ξ} = \frac{(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) (R_{0 {\hat{k}}_{0}}^{S} - 1)}{(β_{0 {\hat{k}}_{0}} - S^{0} σ_{{\hat{k}}_{0}}^{2}) \frac{S^{0}}{b} ((μ_{N} + μ_{{\tilde{k}}_{*}} + γ_{{\tilde{k}}_{*}}) - \frac{δ γ_{{\tilde{k}}_{*}}}{μ_{N} + δ})} = \frac{β_{0 k} S^{0} - \frac{σ_{k}^{2}}{2} {(S^{0})}^{2} - (μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}})}{(S^{0} β_{0 {\hat{k}}_{0}} - S^{02} σ_{{\hat{k}}_{0}}^{2}) D}

(33)

and

\frac{\partial \bar{ξ}}{\partial S^{0}} = \frac{\frac{1}{2} β_{0 k} S^{02} σ_{{\hat{k}}_{0}}^{2} + (μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) (β_{0 {\hat{k}}_{0}} - 2 S^{0} σ_{{\hat{k}}_{0}}^{2})}{{(S^{0} β_{0 {\hat{k}}_{0}} - S^{02} σ_{{\hat{k}}_{0}}^{2})}^{2} D} .

(34)

The numerator of

\frac{\partial \bar{ξ}}{\partial S^{0}}

is

\begin{array}{l} \frac{1}{2} β_{0 k} σ_{{\hat{k}}_{0}}^{2} S^{02} - 2 D σ_{{\hat{k}}_{0}}^{2} S^{0} + (μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) β_{0 {\hat{k}}_{0}} \\ = \frac{1}{2} β_{0 k} σ_{{\hat{k}}_{0}}^{2} {(S^{0} - 2 \frac{(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}})}{β_{0 k}})}^{2} - 2 D {(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}})}^{2} σ_{{\hat{k}}_{0}}^{2} + (μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) β_{0 {\hat{k}}_{0}} . \end{array}

(35)

If

β_{0 {\hat{k}}_{0}} > 2 (μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) σ_{{\hat{k}}_{0}}^{2}

,

\bar{ξ}

is a strictly increasing function of

S^{0}

. Then,

\bar{ξ}

is a decreasing function of

q_{1}

. We can consider a cost function as

C_{1} (\bar{ξ}, q_{1}) = c_{2} \bar{ξ} - c_{0} ln (1 - q_{1}^{2}) - c_{1} ln (1 - q_{1}) .

(36)

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 \in {1, 2, 3}} R_{0 k}^{S} = 0.2941 < 1

and

S^{0} = 2.3810

. It implies

β_{0 k} > S^{0} σ_{k}^{2}

and

\frac{β_{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 \in {1, 2, 3}} R_{0 k}^{S} = 0.2813 < 1

and

{m a x}_{k \in {1, 2, 3}} \frac{β_{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

{\bar{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

\frac{σ_{2}^{2} β_{01}}{σ_{1}^{2} β_{02}} = 1.1200 > 1

,

\frac{σ_{1}^{2}}{(μ_{N} + μ_{1} + γ_{1})} = 0.0104 > 0.0095 = \frac{σ_{2}^{2}}{(μ_{N} + μ_{2} + γ_{2})}

,

\frac{σ_{3}^{2} β_{01}}{σ_{1}^{2} β_{03}} = 0.9333 < 1

, and

\frac{σ_{3}^{2}}{(μ_{N} + μ_{3} + γ_{3})} = 0.0039 > 0.0013 = (\frac{R_{03}^{S} - 1}{R_{01}^{S} - 1}) \frac{σ_{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} \to 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

\bar{ξ} = 3.8743

. If no one is vaccinated at all,

\lim_{q_{1} \to 0} \bar{ξ} = 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

\bar{ξ}

’srelated to different

(q_{0}, q_{1})

combinations are shown inTable 3. All

R_{0}^{S}

’s and

\bar{ξ}

’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 (

\sum_{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

δ = \frac{1}{180}

,

μ_{k} = 0.01 \times \frac{1}{14}

, and

γ_{k} = 0.99 \times \frac{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} = \frac{1}{180}

and

ε_{1} = \frac{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

\frac{0.055}{N}

and

\frac{0.18}{N}

in the table, we assign

(β_{1}, β_{2}, β_{3}) = \frac{1}{1000} \times (0.18, 0.15, 0.06)

and

(σ_{1}, σ_{2}, σ_{3}) = \frac{1}{1000} \times (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) \in Λ

, there is a unique solution

θ (t) \in Λ

for all

t \in [0, τ_{e})

. Let

m_{0} > 0

be a large value such that

θ (0) \in A_{m_{0}}

, where

A_{m_{0}} = {θ (t) \in Λ : \frac{1}{m_{0}} < \min_{1 \leq i \leq 2 K + 4} θ_{i} (t), \max_{1 \leq i \leq 2 K + 4} θ_{i} (t) < m_{0}} .

(A1)

For any

m > m_{0}

, we define the stopping time

τ_{m} = \inf {t \in [0, τ_{e}) : θ (t) \notin A_{m}} .

(A2)

Obviously,

τ_{m}

is increasing with

m

, and

τ_{\infty} = \lim_{m \to \infty} τ_{m} \leq τ_{e}

a.s. If we can prove

τ_{\infty} = \infty

a.s., then

τ_{e} = \infty

a.s. and

θ (t) \in Λ

for all

t \geq 0

a.s.

If

τ_{e} = \infty

a.s. is invalid, there will exist a

T > 0

and an

ϵ > 0

such that

P (τ_{e} \leq T) > ϵ .

(A3)

It means there exists a large enough

m_{1} \geq m_{0}

such that for all

m \geq m_{1}

, we have

P (τ_{m} \leq T) > ϵ .

(A4)

Define a function

W : ℝ_{+}^{K + 4} \to ℝ_{+} \cup {0}

by

W (θ) = (S - 1 - \ln S) + \sum_{k = 1}^{K} (I_{k} - 1 - \ln I_{k}) + (R - 1 - \ln R) + \sum_{i = 1}^{2} (V_{i} - 1 - \ln V_{i}) .

(A5)

By the It

\hat{o}

’s formula,

\begin{matrix} d W & = (1 - \frac{1}{S}) d S + \frac{1}{2 S^{2}} {(d S)}^{2} + \sum_{k = 1}^{K} [(1 - \frac{1}{I_{k}}) d I_{k} + \frac{1}{2 I_{k}^{2}} {(d I_{k})}^{2}] + (1 - \frac{1}{R}) d R \\ + \sum_{i = 1}^{2} [(1 - \frac{1}{V_{i}}) d V_{i} + \frac{1}{2 V_{i}^{2}} {(d V_{i})}^{2}] \\ = L W d t + (1 - \frac{1}{S}) [- S \sum_{k = 1}^{K} σ_{k} I_{k} d B_{k}] + \sum_{k = 1}^{K} (1 - \frac{1}{I_{k}}) [S σ_{k} I_{k} d B_{k}] \\ = L W d t + (\sum_{k = 1}^{K} (I_{k} - S) σ_{k} d B_{k}), \end{matrix}

(A6)

where

\begin{array}{l} L W = (1 - \frac{1}{S}) [b - S \sum_{k = 1}^{K} β_{0 k} I_{k} - (μ_{N} + q_{0}) S + \sum_{i = 1}^{2} ε_{i} V_{i} + δ R] + \frac{1}{2 S^{2}} (S^{2} \sum_{k = 1}^{K} {(σ_{k} I_{k})}^{2}) \\ + \sum_{k = 1}^{K} [(1 - \frac{1}{I_{k}}) (S β_{0 k} I_{k} - (μ_{N} + μ_{k} + γ_{k}) I_{k})] - \sum_{k = 1}^{K} \frac{1}{2 I_{k}^{2}} S^{2} {(σ_{k} I_{k})}^{2} \\ + (1 - \frac{1}{R}) [\sum_{k = 1}^{K} γ_{k} I_{k} - (μ_{N} + δ) R] \\ + (1 - \frac{1}{V_{1}}) [q_{0} S - (μ_{N} + ε_{1} + q_{1}) V_{1}] + (1 - \frac{1}{V_{2}}) [q_{1} V_{1} - (μ_{N} + ε_{2}) V_{2}] \\ \leq [b + \sum_{i = 1}^{2} ε_{i} V_{i} + δ R + \sum_{k = 1}^{K} β_{0 k} I_{k} + (μ_{N} + q_{0})] + \frac{1}{2} \sum_{k = 1}^{K} {(σ_{k}^{2} I_{k})}^{2} \\ + [\sum_{k = 1}^{K} γ_{k} I_{k} + (μ_{N} + δ)] + \sum_{k = 1}^{K} [S β_{0 k} I_{k} + (μ_{N} + μ_{k} + γ_{k})] \\ + \sum_{k = 1}^{K} \frac{1}{2} S^{2} σ_{k}^{2} + [q_{0} S + (μ_{N} + ε_{1} + q_{1})] + [q_{1} V_{1} + (μ_{N} + ε_{2})] \\ \leq b + (K + 4) μ_{N} + δ + q_{0} + q_{1} + ε_{1} + ε_{2} + \sum_{k = 1}^{K} (μ_{k} + γ_{k}) \\ + \frac{b}{μ_{N}} (\sum_{i = 1}^{2} ε_{i} + δ + \sum_{k = 1}^{K} β_{0 k} + \sum_{k = 1}^{K} γ_{k} + q_{0} + q_{1}) + {(\frac{b}{μ_{N}})}^{2} \sum_{k = 1}^{K} (β_{0 k} + σ_{k}^{β 2}) ≔ M, \end{array}

for all

0 < t < T

. Then, we obtain

E W (θ (T \land τ_{m})) \leq W (θ (0)) + M T < \infty .

(A7)

Let

Ψ_{m} = {τ_{m} \leq T}

,

\forall m \geq m_{0}

, and then

P (Ψ_{m}) > ϵ

. For all

ω \in Ψ_{m}

, we have

W (θ (τ_{δ}, ω)) \geq (m - 1 - \ln m) \land (\frac{1}{m} - 1 + \ln m) .

(A8)

It implies

W (θ (0)) + M T \geq E W (θ (T \land τ_{m}) χ_{(ω \in Ψ_{δ})}) \geq ϵ [(m - 1 - \ln m) \land (\frac{1}{m} - 1 + \ln m)],

(A9)

which leads to the contradiction when

m \to \infty

,

\infty > W (θ (0)) + M T \geq \infty .

The proof is complete. □

Appendix B

Proof of Theorem 2.

(a): Using model (12) and Equation (16), we have

${\begin{cases} \frac{S (t) - S (0)}{t} = b - \sum_{k = 1}^{K} β_{0 k} 〈 S I_{k} 〉_{t} - (μ_{N} + q_{0}) {〈 S 〉}_{t} + \sum_{i = 1}^{2} ε_{i} 〈 V_{i} 〉_{t} + δ {〈 R 〉}_{t} \\ - \frac{1}{t} \sum_{k = 1}^{K} σ_{k} \int_{0}^{t} S I_{k} d B_{k} (s), \\ \frac{I_{k} (t) - I_{k} (0)}{t} = β_{0 k} 〈 S I_{k} 〉_{t} - (μ_{N} + μ_{k} + γ_{k}) 〈 I_{k} 〉_{t} + \frac{1}{t} σ_{k} \int_{0}^{t} S I_{k} d B_{k} (s), k = 1, 2, \dots, K, \\ \frac{R (t) - R (0)}{t} = \sum_{k = 1}^{K} γ_{k} 〈 I_{k} 〉_{t} - (μ_{N} + δ) {〈 R 〉}_{t}, \\ \frac{V_{1} (t) - V_{1} (0)}{t} = q_{0} {〈 S 〉}_{t} - (μ_{N} + ε_{1} + q_{1}) 〈 V_{1} 〉_{t}, \\ \frac{V_{2} (t) - V_{2} (0)}{t} = q_{1} 〈 V_{1} 〉_{t} - (μ_{N} + ε_{2}) 〈 V_{2} 〉_{t} . \end{cases}$

(A10)

According to system (A10) of equations, we have

\begin{array}{l} (\frac{S (t) - S (0)}{t} + \sum_{k = 1}^{K} \frac{I_{k} (t) - I_{k} (0)}{t}) + \frac{ε_{1} (μ_{N} + ε_{2}) + ε_{2} q_{1}}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})} \frac{V_{1} (t) - V_{1} (0)}{t} + \frac{ε_{2}}{μ_{N} + ε_{2}} \frac{V_{2} (t) - V_{2} (0)}{t} \\ + \frac{δ}{μ_{N} + δ} \frac{R (t) - R (0)}{t} \\ = b - (μ_{N} + q_{0}) {〈 S 〉}_{t} - \sum_{k = 1}^{K} ((μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}) 〈 I_{k} 〉_{t} + (q_{0} \frac{ε_{1} (μ_{N} + ε_{2}) + ε_{2} q_{1}}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}) {〈 S 〉}_{t} . \end{array}

(A11)

Then,

\begin{matrix} {〈 S 〉}_{t} & = \frac{1}{(μ_{N} + q_{0}) - q_{0} \frac{ε_{1} (μ_{N} + ε_{2}) + ε_{2} q_{1}}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}} (b - \sum_{k = 1}^{K} ((μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}) 〈 I_{k} 〉_{t}) - φ (t) \\ = \frac{b}{μ_{N}} \frac{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2}) + (μ_{N} q_{0} + q_{0} q_{1} + ε_{2} q_{0})} \\ - \frac{1}{μ_{N}} \frac{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2}) + (μ_{N} q_{0} + q_{0} q_{1} + ε_{2} q_{0})} \sum_{k = 1}^{K} ((μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}) 〈 I_{k} 〉_{t} - φ (t) \\ = S^{0} - \frac{S^{0}}{b} \sum_{k = 1}^{K} ((μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}) 〈 I_{k} 〉_{t} - φ (t), \end{matrix}

(A12)

where

\begin{matrix} φ (t) = \frac{S^{0}}{b} [(\frac{S (t) - S (0)}{t} + \sum_{k = 1}^{K} \frac{I_{k} (t) - I_{k} (0)}{t}) & + \frac{ε_{1} (μ_{N} + ε_{2}) + ε_{2} q_{1}}{(μ_{N} + ε_{1} + q_{1}) (μ_{N} + ε_{2})} \frac{V_{1} (t) - V_{1} (0)}{t} \\ + \frac{ε_{2}}{μ_{N} + ε_{2}} \frac{V_{2} (t) - V_{2} (0)}{t} + \frac{δ}{μ_{N} + δ} \frac{R (t) - R (0)}{t}] . \end{matrix}

(A13)

Note that

\lim_{t \to \infty} \frac{N (t)}{t} \leq \lim_{t \to \infty} \frac{1}{t} \frac{b}{μ_{N}} = 0

, we have

\lim_{t \to \infty} φ (t) = 0, a . s .

(A14)

Applying It

\hat{o}

’s formula to the equations of

I_{k}

in model (12) and using (16) again, we obtain

\frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} = β_{0 k} {〈 S 〉}_{t} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} 〈 S^{2} 〉_{t} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s), a . s .

(A15)

By the Cauchy–Schwarz inequality,

〈 S^{2} 〉_{t} \geq {〈 S 〉}_{t}^{2} a . s .

Then,

\frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} \leq β_{0 k} {〈 S 〉}_{t} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {〈 S 〉}_{t}^{2} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s), a . s .

(A16)

Now, using (A12) and (A16),

\begin{array}{l} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} \leq β_{0 k} [S^{0} - \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - φ (t)] \\ - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} [{(S^{0})}^{2} + {(\frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t})}^{2} + φ^{2} (t) \\ - 2 \frac{S^{0}}{b} (S^{0} - φ (t)) \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - 2 S^{0} φ (t)] \\ + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s) \end{array}

(A17)

We define the last term of Equation (A17) as

M_{k} (t) = σ_{k} \int_{0}^{t} S d B_{k} (s) .

(A18)

Then,

〈 M_{k}, M_{k} 〉_{t} = \int_{0}^{t} σ_{k}^{2} S^{2} (r) d r \leq σ_{k}^{2} {(\frac{b}{μ_{N}})}^{2} t

, and we have

\underset{t \to \infty}{limsup} \frac{〈 M_{k}, M_{k} 〉_{t}}{t} \leq σ_{k}^{2} {(\frac{b}{μ_{N}})}^{2} < \infty .

(A19)

By Lemma 1,

\lim_{t \to \infty} \frac{M_{k} (t)}{t} = 0, a . s .

(A20)

Plugging

M_{k} (t)

into Equation (A17),

\begin{array}{l} \frac{\ln I_{k} (t)}{t} \leq \frac{\ln I_{k} (0)}{t} + β_{0 k} [S^{0} - \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - φ (t)] \\ - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} [S^{02} + {(\frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t})}^{2} + φ^{2} (t) \\ - 2 \frac{S^{0}}{b} (S^{0} - φ (t)) \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - 2 S^{0} φ (t)] + \frac{M_{k} (t)}{t} \\ \leq \frac{\ln I_{k} (0)}{t} + β_{0 k} [S^{0} - φ (t)] - (μ_{N} + μ_{k} + γ_{k}) \\ - \frac{σ_{k}^{2}}{2} [S^{02} + {(\frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t})}^{2} + φ^{2} (t) \\ + 2 φ (t) \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - 2 S^{0} φ (t)] \end{array} - (β_{0 k} - S^{0} σ_{k}^{2}) \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} + \frac{M_{k} (t)}{t} .

(A21)

By the assumption

β_{0 k} > S^{0} σ_{k}^{2}

, we have

\begin{matrix} \frac{\ln I_{k} (t)}{t} \leq \frac{\ln I_{k} (0)}{t} & + β_{0 k} S^{0} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} S^{02} - β_{0 k} φ (t) \\ - \frac{σ_{k}^{2}}{2} [2 φ (t) \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - 2 S^{0} φ (t)] . \end{matrix}

(A22)

As

R_{0 k}^{S} - 1 = R_{0 k} - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2} - 1 < 0

, we have

\begin{matrix} \underset{t \to \infty}{limsup} \frac{\ln I_{k} (t)}{t} & \leq β_{0 k} S^{0} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {(S^{0})}^{2} \\ = (μ_{N} + μ_{k} + γ_{k}) [R_{0 k} - 1 - \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} {(S^{0})}^{2}] < 0, a . s . \end{matrix}

(A23)

(b): If $\frac{β_{0 k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} < σ_{k}^{2}$ , by Equation (A16),

$\begin{matrix} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} & \leq β_{0 k} {〈 S 〉}_{t} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {〈 S 〉}_{t}^{2} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s) \\ = - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {({〈 S 〉}_{t} - \frac{β_{0 k}}{σ_{k}^{2}})}^{2} + \frac{β_{0 k}^{2}}{2 σ_{k}^{2}} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s) \\ \leq - (μ_{N} + μ_{k} + γ_{k}) + \frac{β_{0 k}^{2}}{2 σ_{k}^{2}} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s) . \end{matrix}$

(A24)

It implies

\underset{t \to \infty}{limsup} \frac{\ln I_{k} (t)}{t} \leq - (μ_{N} + μ_{k} + γ_{k}) + \frac{β_{0 k}^{2}}{2 σ_{k}^{2}} < 0, a . s .

(A25)

The proof is complete. □

Appendix C

Proof of Theorem 4.

(a): Using (A12) and (A16),

$\begin{array}{l} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} = β_{0 k} {〈 S 〉}_{t} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} 〈 S^{2} 〉_{t} + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s) \\ \geq β_{0 k} S^{0} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {(\frac{b}{μ_{N}})}^{2} - β_{0 k} \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \end{array} + β_{0 k} φ (t) + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s),$

(A26)

We have

β_{0 k} \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \geq - \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} + β_{0 k} S^{0} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {(\frac{b}{μ_{N}})}^{2} + β_{0 k} φ (t) + \frac{1}{t} σ_{k} \int_{0}^{t} S d B_{k} (s),

(A27)

k = 1, 2, \dots, N

. An upper bound of the left side of inequality (A27) is

β_{0 k} \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) \frac{b}{μ_{N}},

implying

\underset{t \to \infty}{limsup} (- \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t})

has an upper bound. As

I_{k} (t) \in (0, \frac{b}{μ_{N}})

, we have

\underset{t \to \infty}{liminf} (- \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t}) \geq 0, a . s .

(A28)

Then,

\underset{t \to \infty}{liminf} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \geq \frac{(μ_{N} + μ_{k} + γ_{k}) (R_{0 k}^{*} - 1)}{β_{0 k} \frac{S^{0}}{b}},

(A29)

k = 1, 2, \dots, N

. We set

{\hat{k}}^{*} = \arg \max_{k : R_{0 k}^{*} > 1} {\frac{(μ_{N} + μ_{k} + γ_{k}) (R_{0 k}^{*} - 1)}{β_{0 k}}}

. Then

\underset{t \to \infty}{liminf} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \geq \frac{(μ_{N} + μ_{{\hat{k}}^{*}} + γ_{{\hat{k}}^{*}}) (R_{0 {\hat{k}}^{*}}^{*} - 1)}{β_{0 {\hat{k}}^{*}} \frac{S^{0}}{b}}, a . s .

(A30)

Setting

{\tilde{k}}^{*} = \arg \max_{k : R_{0 k}^{*} > 1} {(μ_{N} + μ_{k} + γ_{k}) - \frac{δ γ_{k}}{μ_{N} + δ}}

, Equation (22) is held.

(b): Using (A12) and (A17),

$\begin{matrix} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} \leq β_{0 k} S^{0} & - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {(S^{0})}^{2} + ϕ (t) \\ - (β_{0 k} - S^{0} σ_{k}^{2}) \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} + \frac{M_{k} (t)}{t}, \end{matrix}$

(A31)

where

$\begin{matrix} ϕ (t) = - β_{0 k} φ (t) \\ - \frac{σ_{k}^{2}}{2} [φ^{2} (t) + 2 φ (t) \frac{S^{0}}{b} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} - 2 S^{0} φ (t)] . \end{matrix}$

(A32)

By Equation (A14),

\lim_{t \to \infty} φ (t) = 0

implies

\lim_{t \to \infty} ϕ (t) = 0, a . s .

Then, we have

\begin{matrix} (β_{0 k} - S^{0} σ_{k}^{2}) \frac{S^{0}}{b} & \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \\ \leq - \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} + β_{0 k} S^{0} - (μ_{N} + μ_{k} + γ_{k}) - \frac{σ_{k}^{2}}{2} {(S^{0})}^{2} + ϕ (t) + \frac{M_{k} (t)}{t} . \end{matrix}

(A33)

By assumption, there exists some

k_{0} \in {1, 2, \dots, K}

such that

\underset{t \to \infty}{liminf} 〈 I_{k_{0}} 〉_{t} > 0

,

β_{0 k_{0}} > S^{0} σ_{k_{0}}^{2}

, and

R_{0 k_{0}}^{S} \geq R_{0 k_{0}}^{*} > 1

. With

I_{k_{0}} \leq \frac{b}{μ_{N}}

, we obtain

\lim_{t \to \infty} \frac{\ln I_{k_{0}} (t) - \ln I_{k_{0}} (0)}{t} = 0, a . s .

Thus,

\begin{matrix} (β_{0 k_{0}} - S^{0} σ_{k_{0}}^{2}) \frac{S^{0}}{b} \underset{t \to \infty}{limsup} \sum_{k^{'} = 1}^{K} ((μ_{N} & + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \\ \leq β_{0 k_{0}} S^{0} - (μ_{N} + μ_{k_{0}} + γ_{k_{0}}) - \frac{σ_{k}^{2}}{2} {(S^{0})}^{2} . \end{matrix}

(A34)

We have

\underset{t \to \infty}{limsup} \sum_{k^{'} = 1}^{K} ((μ_{N} + μ_{k^{'}} + γ_{k^{'}}) - \frac{δ γ_{k^{'}}}{μ_{N} + δ}) 〈 I_{k^{'}} 〉_{t} \leq \frac{(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) (R_{0 {\hat{k}}_{0}}^{S} - 1)}{(β_{0 {\hat{k}}_{0}} - S^{0} σ_{{\hat{k}}_{0}}^{2}) \frac{S^{0}}{b}},

(A35)

where

{\hat{k}}_{0} : = {argmax}_{k_{0} : R_{0 k_{0}}^{S} > 1} {\frac{(μ_{N} + μ_{k_{0}} + γ_{k_{0}}) (R_{0 k_{0}}^{S} - 1)}{(β_{0 k_{0}} - S^{0} σ_{k_{0}}^{2})}}

.

By the definition of

{\tilde{k}}_{*}

, we have

\underset{t \to \infty}{limsup} ((μ_{N} + μ_{{\tilde{k}}_{*}} + γ_{{\tilde{k}}_{*}}) - \frac{δ γ_{{\tilde{k}}_{*}}}{μ_{N} + δ}) \sum_{k^{'} = 1}^{K} 〈 I_{k^{'}} 〉_{t} \leq \frac{(μ_{N} + μ_{{\hat{k}}_{0}} + γ_{{\hat{k}}_{0}}) (R_{0 {\hat{k}}_{0}}^{S} - 1)}{(β_{0 {\hat{k}}_{0}} - S^{0} σ_{{\hat{k}}_{0}}^{2}) \frac{S^{0}}{b}}

(A36)

and then Equation (23) is held.

The proof is complete. □

Appendix D

Proof of Theorem 5.

Using Equation (A15),

\begin{array}{l} β_{0 {\bar{k}}^{*}} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} - β_{0 k} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} \\ = (β_{0 {\bar{k}}^{*}} β_{0 k} - β_{0 k} β_{0 {\bar{k}}^{*}}) {〈 S 〉}_{t} - β_{0 {\bar{k}}^{*}} (μ_{N} + μ_{k} + γ_{k}) + β_{0 k} (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}) \\ - (β_{0 {\bar{k}}^{*}} \frac{σ_{k}^{2}}{2} - β_{0 k} \frac{σ_{{\bar{k}}^{*}}^{2}}{2}) 〈 S^{2} 〉_{t} + \frac{1}{t} β_{0 {\bar{k}}^{*}} σ_{k} \int_{0}^{t} S d B_{k} (s) - \frac{1}{t} β_{0 k} σ_{{\bar{k}}^{*}} \int_{0}^{t} S d B_{{\bar{k}}^{*}} (s) \\ = (μ_{N} + μ_{k} + γ_{k}) (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}) \times \\ [\frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} - \frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} - (\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{2 (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}) 〈 S^{2} 〉_{t}] \end{array} + \frac{1}{t} β_{0 {\bar{k}}^{*}} σ_{k} \int_{0}^{t} S d B_{k} (s) - \frac{1}{t} β_{0 k} σ_{{\bar{k}}^{*}} \int_{0}^{t} S d B_{{\bar{k}}^{*}} (s)

(A37)

By assumptions,

\begin{matrix} \lim_{t \to \infty} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} = \lim_{t \to \infty} (β_{0 {\bar{k}}^{*}} {〈 S 〉}_{t} - (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}) - \frac{σ_{{\bar{k}}^{*}}^{2}}{2} 〈 S^{2} 〉_{t} + \frac{1}{t} σ_{{\bar{k}}^{*}} \int_{0}^{t} S d B_{{\bar{k}}^{*}}) \\ = (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}) \lim_{t \to \infty} (\frac{β_{0 {\bar{k}}^{*}} {〈 S 〉}_{t}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} - \frac{σ_{{\bar{k}}^{*}}^{2}}{2 (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} 〈 S^{2} 〉_{t} - 1) = 0, a . s . \end{matrix}

(A38)

There exists a large enough

T_{0} > 0

such that for all

t \geq T_{0}

, and all

ω \in Ω_{0} \subseteq Ω

with

P (Ω_{0}) = 1

, we have

〈 S^{2} 〉_{t} (ω) = \frac{2 β_{0 {\bar{k}}^{*}} {〈 S 〉}_{t} (ω)}{σ_{{\bar{k}}^{*}}^{2}} - \frac{2 (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}{σ_{{\bar{k}}^{*}}^{2}} + ϵ_{t} (ω),

(A39)

where

| ϵ_{t} (ω) | \geq 0

and

| ϵ_{t} (ω) | \to 0

when

t \to \infty

.

Substituting (A39) to the first term on the right-hand side of (A37), we have

\begin{array}{l} β_{0 {\bar{k}}^{*}} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} - β_{0 k} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} \\ = β_{0 {\bar{k}}^{*}} \frac{(μ_{N} + μ_{k} + γ_{k}) (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}{σ_{{\bar{k}}^{*}}^{2}} [- \frac{σ_{{\bar{k}}^{*}}^{2}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \\ - (\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}) {〈 S 〉}_{t} (ω) + \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})}] + ϵ_{t}^{*} (ω), \end{array}

(A40)

where

\begin{array}{l} ϵ_{t}^{*} (ω) = ϵ_{t} (ω) (\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{2 (μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{2 (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}) \\ + \frac{1}{t} β_{0 {\bar{k}}^{*}} σ_{k} \int_{0}^{t} S d B_{k} (s, ω) - \frac{1}{t} β_{0 k} σ_{{\bar{k}}^{*}} \int_{0}^{t} S d B_{{\bar{k}}^{*}} (s, ω) . \end{array}

(A41)

Now, if

\frac{β_{0 {\bar{k}}^{*}}}{σ_{{\bar{k}}^{*}}^{2}} \geq \frac{β_{0 k}}{σ_{k}^{2}}

, then

\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} > 0 .

As

S_{t} \geq 0 a . s .

and

\frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} > \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})}

, the right hand side of Equation (A40) is negative when

ϵ_{t}^{*} (ω)

is dropped. Then,

\underset{t \to \infty}{limsup} β_{0 {\bar{k}}^{*}} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} - β_{0 k} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} < 0 a . s .

(A42)

For

\lim_{t \to \infty} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} = 0, a . s .

, we have

β_{0 {\bar{k}}^{*}} \underset{t \to \infty}{limsup} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} < 0, a . s .

If

\frac{β_{0 {\bar{k}}^{*}}}{σ_{{\bar{k}}^{*}}^{2}} < \frac{β_{0 k}}{σ_{k}^{2}}

, then

\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} < 0 .

By Equation (A12), there exists a

T_{1} > 0

large enough such that

{〈 S 〉}_{t} \geq S^{0}, a . s .

\forall t \geq T_{1}

. Thus, we have

Ω_{1} \subseteq Ω

,

P (Ω_{1}) = 1

, such that for all

ω \in Ω_{0} \cap Ω_{1}

and

t > \max {T_{0}, T_{1}}

,

\begin{array}{l} β_{0 {\bar{k}}^{*}} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} - β_{0 k} \frac{\ln I_{{\bar{k}}^{*}} (t) - \ln I_{{\bar{k}}^{*}} (0)}{t} \\ \leq β_{0 {\bar{k}}^{*}} \frac{(μ_{N} + μ_{k} + γ_{k}) (μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}{σ_{{\bar{k}}^{*}}^{2}} [- \frac{σ_{{\bar{k}}^{*}}^{2}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \\ - (\frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} - \frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}) S^{0} + \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})}] + ϵ_{t}^{*} (ω) . \end{array}

(A43)

The part in the rectangle bracket on the right-hand side of Equation (A43) is

\begin{array}{l} - \frac{σ_{{\bar{k}}^{*}}^{2}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} + (\frac{β_{0 k}}{μ_{N} + μ_{k} + γ_{k}} \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} - \frac{β_{0 {\bar{k}}^{*}}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})}) S^{0} - \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} \\ = (\frac{β_{0 k} S^{0}}{μ_{N} + μ_{k} + γ_{k}} - 1) \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} - (\frac{β_{0 {\bar{k}}^{*}} S^{0}}{μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}}} - 1) \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} \\ = (R_{0 k}^{S} - 1) \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})} - (R_{0 {\bar{k}}^{*}}^{S} - 1) \frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} . \end{array}

(A44)

Then, if

\frac{σ_{k}^{2}}{(μ_{N} + μ_{k} + γ_{k})} > (\frac{R_{0 k}^{S} - 1}{R_{0 {\bar{k}}^{*}}^{S} - 1}) \frac{σ_{{\bar{k}}^{*}}^{2}}{(μ_{N} + μ_{{\bar{k}}^{*}} + γ_{{\bar{k}}^{*}})}

, the right hand side of Equation (A43) is negative when

ϵ_{t}^{*} (ω)

is dropped.

Again, it implies

\underset{t \to \infty}{limsup} \frac{\ln I_{k} (t) - \ln I_{k} (0)}{t} < 0, a . s .

The proof is complete. □

References

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]
Anderson, R.M.; May, R.M. Population biology of infectious diseases: Part I. Nature 1979, 280, 361–367. [Google Scholar] [CrossRef] [PubMed]
May, R.M.; Anderson, R.M. Population biology of infectious diseases: Part II. Nature 1979, 280, 455–461. [Google Scholar] [CrossRef] [PubMed]
Capasso, V.; Serio, G. A generalization of the Kermack-McKendrick deterministic epidemic model. Math. Biosci. 1978, 42, 43–61. [Google Scholar] [CrossRef]
Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: Boston, MA, USA, 2015. [Google Scholar]
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]
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]
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]
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]
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]
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]
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]
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]
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]
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]
Li, J.; Ma, Z. Global analysis of SIS epidemic models with variable total population size. Math. Comput. Model. 2004, 39, 1231–1242. [Google Scholar]
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]
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]
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]
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).
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]
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]
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]
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]
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]
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]
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]
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]
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 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}

.

Figure 3. The simulated numbers of infected individuals under condition (a) of Theorem 2.

Figure 4. The simulated numbers of infected individuals under condition (b) of Theorem 2.

Figure 5. The simulated numbers of infected individuals under the conditions of the disease persistence.

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.

Figure 7. The simulated numbers of infected individuals under the conditions of Table 1 in [22].

Table 1. Description of variables and parameters in model (1).

Variable/Parameter	Description
$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.

Parameter	PRCC	Parameter	PRCC
$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

\bar{ξ}

’s.

Table 3. The different

(q_{0}, q_{1})

combinations and related

R_{0}^{S}

’s and

\bar{ξ}

’s.

$(q_{0}, q_{1})$	$R_{0}^{S}$	$\bar{ξ}$
(0.1, 0.1)	8.6662	5.7327
(0.1, 0.5)	8.1300	5.6695
(0.1, 0.9)	7.9410	5.6458
(0.5, 0.1)	4.5833	4.9769
(0.5, 0.5)	3.9022	4.7210
(0.5, 0.9)	3.6912	4.6240
(0.9, 0.1)	3.1151	4.2961
(0.9, 0.5)	2.5669	3.8537
(0.9, 0.9)	2.4041	3.6848

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

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 Menu

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

Abstract

1. Introduction

2. The Model

2.1. The Equilibrium Point and the Basic Reproduction Number

2.2. Sensitivity Analysis of the Basic Reproduction Number to Parameters

2.3. The Stochastic Model

3. The Asymptotic Behavior of the Stochastic Model

3.1. Existence of Unique Positive Solution

3.2. Extinction

3.3. Persistence

3.4. Vaccination Rates

4. Numerical Examples

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Appendix A

Appendix B

Appendix C

Appendix D

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI