The Global dynamics of a SIR model considering competitions among multiple strains in patchy environments

1.   Introduction

Differential equation models, especially ordinary differential equation models, as mathematical tools for studying the laws of biological phenomena over time, have been widely used in the study of the spread of infectious diseases ^[1], the growth of biological populations ^[2], and the laws of gene expression ^[3]. Among them, in the study of infectious diseases, the ordinary differential equation model can be established by classifying the susceptible, the infected and the recovered into different variables ^[1]. By analyzing the dynamics of the model and constructing the basic reproductive number of the model, the persistence and extinction conditions of infectious diseases can be obtained. Meanwhile, based on actual data, the infectious disease ordinary differential equation model can predict the future trend of the disease, and then provide a corresponding basis for the formulation of epidemic prevention policies ^[4].

In the study of infectious disease models, a more special category is the multi-strain model. The research object of this type of model is the spread and competition between different strains during the spread of infectious diseases. For infectious diseases caused by viruses, such as influenza and Covid-19, this model often shows that different viruses are simultaneously spreading among people due to the mutation of the pathogen. Bichara et al. ^[5] construct a series of Lyapunov functions to analyze the global dynamics of the SIS, SIR and MSIR models with multi-strains, and give the conditions whether the disease-free equilibrium could be globally stable or not. Moreover, Meehan et al. ^[6] consider a multiple strain model with mutation, and numerically find the conditions that could influence the coexistence and competition exclusion of that model. Li et al. ^[7] use a multi-strain partial differential equation model to analyze competitions of different strains of the virus under different drug environments, and find the relationship between the basic reproductive number and the drug parameters. All of these works show the importance on analyzing multiple strain models.

Another common type of infectious disease model is the patch model ^[8,9]. The model assumes that, during the spread of infectious diseases, populations in the same patch can infect each other, and populations among different patches can only migrate but not infect each other. This model is suitable for infectious diseases spreading globally, and different countries and regions can be regarded as a patch. In the study of this type of model, Wang et al. ^[10] construct a patchy model to analyze the population dispersal and its effect on the transmission of disease. While, Zhang et al. ^[11] consider a more general patchy model and analyze its global behaviors with the properties of a monotone system. Li et al. ^[12] obtain the global stability of an autonomous system with patches by constructing a Lyapunov function. All of these works analyze the dynamics of different patch models by introducing the basic reproductive number which is the main threshold value to determine whether the epidemic could persist or not.

Some works consider multiple strains and patchy environments. In reference ^[13], Arino et al. consider a SEIR model with multiple strains and patches and use simulation method to analyze the dynamic behaviors of that model. Marva et al. ^[14] use the SIS model without birth rate to analyze the multiple strains and patches effect on the epidemic. In the analysis of that model, the authors consider a special condition that the model could be written as a slow-fast system. Qiu et al. ^[15] consider a multiple strains model that the susceptible population could be neglected in patches and find the persistence and coexistence conditions of that system. In this work, we construct a novel SIR model considering both patches and strains and analyze its global dynamic without any additional assumptions. The model is as follows.

Let the variables $S_i(t)$ illustrate the number of susceptible population at the patch $i$, and $I_{i, k}(t)$ note the number of infected population of the $k$ strain at the patch $i$. Then we have the equation of multiple strains compete in patchy environment.

where $k = 1, 2, ..., N_S$, and $N_S$ is the strain number. Meanwhile, $i = 1, 2, .., N_P$, and $N_P$ is the total number of patch. $m_{ij}$ is the migration rate from the path $i$ to $j$, and the parameters $m_{i} = \sum_{j\neq i} m_{ij}$. $0 < \alpha < 1$ is the rate that describe the ability the illness could slow down the migration rate. $\lambda_i$ is the constant population growth rate, and $d_i$ is the natural mortality at the patch $i$. $\beta_k$, $\mu_k$ and $\gamma_k$ are the infection rate, disease fatality rate and recovery rate of the $k$ strain. All of the parameters are positive. In the next two sections we will prove the basic properties and analyze the global dynamics of this model.

2.   Basic properties

The right-hand side of the sub-system about the susceptible population of system (1.1) could be written in the matrix form

where $S(t) = (S_1(t), S_2(t), \cdots, S_{N_P}(t))^T$, and $I_k(t) = (I_{1, k}(t), I_{2, k}(t), ..., I_{N_P, k}(t))^T$ are column vectors. The matrix,

is the migration matrix. $\Lambda(I_k(t))$ is the diagonal matrix with the diagonal elements $(I_{1, k}(t), I_{2, k}(t), ..., I_{N_P, k}(t))$, $L = (\lambda_1, \lambda_2, \cdots, \lambda_{N_P})$ is the column vectors. $K$ is the diagonal matrix with the diagonal elements $(d_1, d_2, \cdots, d_{N_P})$. As $m_i = \sum_{j\neq i} m_{ij}$ and $m_{ij} > 0$, the matrix $M^T$ is the diagonally dominant matrix, and $M^T-K$ is the strictly diagonally dominant matrix. Thus, by Levy–Desplanques Theorem ^[16], $M^T-K$ is invertible. As the diagonal elements of $M^T-K$ are less than $0$, by Gershgorin circle Theorem ^[16], the real part of eigenvalues of $M^T-K$ are less than $0$. Therefore, $M-K$ is the invertible matrix with the real part of its eigenvalues less than $0$ (the matrix and its transpose share the same eigenvalues). Hence, the following lemma could be given.

Lemma 1. If $M^T$ is a diagonally dominant matrix with the diagonal elements smaller than $0$, and $K$ is a diagonal matrix with the diagonal elements greater than $0$. Then the matrix $M-K$ is invertible and the real parts of its eigenvalues are less than $0$.

The system without the disease could be given as

Thus, this system has the only equilibrium $S^* = (K-M)^{-1}L$ (Lemma 1). As the diagonal elements of the diagonally dominant matrix, $K-M^T$ are larger than $0$, while the non-diagonal elements of that matrix are less than $0$, which implies the matrix $K-M^T$ is an M-matrix ^[17]. By the inverse-positive property of the M-matrix ^[17], the inverse matrix $(K-M^T)^{-1}$ is positive. Hence, the matrix $(K-M)^{-1}$ is positive, and as the vector $L$ is positive, the equilibrium $S^*$ is also positive. Thus, the disease-free equilibrium $E_0$ could be given inn a matrix form as $(S^*, O, .., O, O)$, where $O$ is the $N_P\times 1$ zero vector and the number of $O$ is $N_S+1$.

Similarly, $I_k$ also could be written as the matrix form,

where the matrix $D_k$ is the diagonal matrix with the diagonal elements as $(\mu_k+ \gamma_k + d_1, \cdots, \mu_k+ \gamma_k + d_{N_P}).$ $\Lambda(S(t))$ is the diagonal matrix with the diagonal elements as $(S_1(t), S_2(t), \cdots, S_{N_P}(t))$. Follows the definitions of the basic reproductive number in patchy system ^[11], the definitions of the basic reproductive number of the strain $k$, could be given as radius

where the function $r(\cdot)$ denote the spectral radius of the matrix.

The positivity and boundedness of the system (1.1) could be given as the following lemma.

Lemma 2 (Boundedness).Solutions of model (1.1) with positive initial conditions are positive and bounded.

Proof. The positive could be proved by contradiction. Assuming the solution of positive initial conditions could go to negative, then there exist a time for $S_i$, $I_{i, k}$ or $R_i$, $i = 1, 2, .., N_P$, $k = 1, 2, \cdots, N_S$ become zero. When $t_0$ is the first time that $S_i = 0$ and other variables greater than $0$, then we have

which follows, there exist a small constant $\delta > 0$, such that $S_i(t) < 0$, when $t\in(t_0-\delta, t_0)$, which leads to the contradiction. Thus, $S_i$ couldn't be the first time for the system touch the boundary of $R^{n+}$. Similarly, we have,

and,

Thus, none of the variables of the system (1.1) could across the boundary of $R^{n+}$. Thus we proved the positivity of the system (1.1).

Then, we will prove the boundedness of the system (1.1). Assuming $N_i(t) = S_i(t) + \sum_k I_{i, k}(t) + R_i(t)$, $i = 1, 2, \cdots, N_P$, then we have $N_i(t) > 0$ and

which is a linear system and right hand side of that system satisfy the Lipschitz condition. Thus the solution of this system is in $C^1$ (page 148 of reference ^[18]), which follows the solution couldn't be infinity at finite time. The Lemma 1 implies the linear system (2.1) have the unique positive equilibrium $S^*$ and the real part of the eigenvalues with respect to that equilibrium are less than $0$, which implies the global stability of that equilibrium, thus the solution of system (2.1) couldn't be infinity when the time goes to infinity. Therefore, all of the solutions of the linear system (2.1) with positive initial condition are bounded. Thus, the system (1.1) is also bounded.

Remark: We can also get the invariant set of the system (1.1) by the properties of cooperating system. The parameters $m_{ji} > 0$ implies the system (2.1) is a cooperating system (^[19], chapter 3, remark 1.1). The monotonicity of the cooperation system (^[19], chapter 3, proposition 2.1) and the global stability of the equilibrium $S^*$ implies that all of the trajectories that comes from a vector $v\in A$, $A = \{(v_1, v_2, \cdots, v_{N_p})|0 < v_i < (S^*)_i, i = 1, 2, \cdots, N_P\}$, will stay in $A$. Thus, the set $A_I = \{(S, I_1, I_2, \cdots, I_{N_S}, R)|S\in A, I_i\in A, i = 1, 2, \cdots, N_S, R\in A\}$, is the invariant set of the solutions of system (1.1). In the rest of this manuscript, we only consider the solutions in the invariant set $A_I$.

3.   Global dynamics

In this section, we mainly talk about the global dynamics of system (1.1). And in the following proofs, the comparison symbols($<$ and $>$) between two vectors compare all of the corresponding elements of the vectors. Which means $(a_1, a_2, \cdots, a_n) < (b_1, b_2, \cdots, b_n)$ is $a_i < b_i, i = 1, 2, \cdots, n$, and same for the signs $>$, $\le$ and $\ge$.

Theorem 1. For model (1.1), if $R_{0, k} < 1$, $k = 1, 2, \cdots, N_S$, then the disease-free equilibrium $E_0$ is globally asymptotically stable.

Proof. Take the derivative of $\frac{dI_{i, k}(t)}{dt}$, with respect to $I_{j, k}$, $i\neq j$, we could have

Hence, the sub-system of system (1.1) that constructed with $I_{i, k}, i = 1, 2, \cdots, N_P$ is a cooperation system (^[19], chapter 3, remark 1.1).

As we only consider the solution in $A_I$ and $S(t) < S^*$, thus,

right hand side of the above inequation is the linear system $I'_k = (\alpha M + \beta_k\Lambda(S^*) -D_k)I_k$, which have the unique equilibrium $(0, 0, \cdots, 0)^T$ ($N_P$ zeros), and the condition $R_{0, k} = r\Big(exp(\alpha M + \beta_k\Lambda(S^*)-D_k)\Big) < 1$ implies the global stability of that zeros equilibrium (Similar to the proof in Lemma 2). Thus, by the comparison of the cooperative system (comes from the monotonicity of cooperative system and the comparison Theorem in the page 112 of reference ^[18]), we could have $\lim_{t\to +\infty}I_k(t) = (0, 0, \cdots, 0)^T$. Thus, $\lim_{t\to +\infty}I_{i, k}(t) = 0$, $i = 1, 2, \cdots, N_P$, $k = 1, 2, \cdots, N_S$. Which follows, $\forall \delta > 0$, there exist a time $T_1 > 0$, such that $\forall t > T_1$, we have $|\beta_{k}I_{i, k}(t)| < \delta$. Thus, $\forall \delta > 0$ and $t > T_1$, the following equations hold,

Choosing a constant $\delta > 0$ and considering the linear system with the variables as $X(t) = (x_1(t), x_2(t), \cdots, x_{N_P})^T$,

which could be written as the matrix form $X' = L+(M-K-\delta E)X$. The Lemma 1 implies the system (3.1) have the unique equilibrium $(K+\delta E-M)^{-1}L$, and the real part of eigenvalues of $M-K-\delta E$ is less than $0$. Which leads to the global stability of the equilibrium $(K+\delta E-M)^{-1}L$ of system (3.1). Thus, $\forall \epsilon > 0$, there exist a time $T_2$, when $t > T_2$, the solution $X(t)$ of the system (3.1) satisfies $|X(t)-(K+\delta E-M)^{-1}L| < \epsilon v_E$, where $v_E$ is the vector with all of its elements are $1$.

Take the derivative of $\frac{dS_i(t)}{dt}$, with respect to $S_j$, $i\neq j$, we could have

Hence, the sub-system of system (1.1) that constructed with $S_i, i = 1, 2, \cdots, N_P$ is a cooperation system. Thus, by the comparison of cooperative system, $\forall\epsilon > 0$ and $\delta > 0$, when $T > \max(T_1, T_2)$ the solution $S(t) = (S_1(t), S_2(t), \cdots, S_{N_p}(t))$ of system (1.1) satisfy the condition, $S(t)\ge (K+\delta E-M)^{-1}L-\epsilon v_E$. From the arbitrary of $\delta$ and $\epsilon$, and $S(t)\le S^* = (K-M)^{-1}L$ ($S(t)\in A$), we have $\lim_{t\to+\infty}S(t) = S^*$. Thus, we proved the global attractivity of the system (1.1) with the basic reproductive numbers smaller than 1.

Then, we go to prove the local stability. By calculating the Jacobin matrix of the disease-free equilibrium of system (1.1), we could find the eigenvalues of that Jacobin matrix calculated from $N_S+2$ sub-matrices. These sub-matrices could be listed as,

Which calculated from the susceptible population, infected population and recovered population. From the conditions $R_{0, k} < 1$, we have the eigenvalues of $\alpha M +\beta_1\Lambda(S^*) -D_1$ is smaller than $0$ (By the definition of the exponential of matrix, it is easy to check that if $\mu$ is the eigenvalue of matrix $M$, then $e^{\mu}$ is the eigenvalue of $\exp(M)$). By Lemma 1, the eigenvalues of $M-K$ are also less than $0$. Thus the local stability of $E_0$ holds. Hence, we complete the proof.

Remark: From the above proof, it is easy to see when one of the real part of the eigenvalues of $\alpha M +\beta_k\Lambda(S^*) -D_k$, $k = 1, 2, \cdots, N_S$ is larger than $0$, the disease free equilibrium $E_0$ will become unstable. It also can be deduced that $E_0$ will be unstable if there is any one $R_{0, k}, k = 1, 2, ···, N_S$ larger than $1$.

Lemma 3. If $M = (m_{ij})$ is a $n\times n$ real matrix, and $\Lambda$ is the diagonal $n\times n$ matrix satisfying $\Lambda = \lambda E$. Then the matrix $\exp(M)$ and $\exp(M+\Lambda)$ share the same eigenvectors, and the corresponding eigenvalues of that two matricessatisfy $\mu = e^{\lambda}\mu_M$, where $\mu_M$ is the eigenvalue of $\exp(M)$ and$\mu$ is the corresponding eigenvalue of $\exp(M+\Lambda)$.

Proof. As $M\Lambda = \Lambda M = \lambda M$, we have $\exp(M+\Lambda) = \exp(M)\exp(\Lambda) = e^{\lambda}\exp(M)$. Let $v$ be the eigenvector of eigenvalue $\mu_M$ of matrix $\exp(M)$, then we have $\exp(M+\Lambda)v = e^{\lambda}\exp(M)v = e^{\lambda}\mu_Mv$, which implies $\mu = e^{\lambda}\mu_M$.

Theorem 2 (Persistence).For model (1.1), if $R_{0, k} > 1$, $k = 1, 2, \cdots, N_S$, then thereexists a positive constant $\delta^* > 0$, such that for at least onestrain $k$ at patch $i$ satisfies $\limsup_{t\to+\infty} I_{i, k}(t)\ge\delta^*$.

Proof. We use the contradictory method to prove this theorem. Assume $\forall \epsilon > 0$, there exits a time $T > 0$, such that $I_{i, k}(t) < \epsilon$. Then similar as the global attractivity proofs in Theorem 1, we could have, $\forall \epsilon > 0$, there exits a time $T > 0$, such that $|S(t)-S^*| < \epsilon v_E$.

As $R_{0, k} > 1$, then the inequality $r\Big(\exp(\alpha M + \beta_k\Lambda(S^*)-D_k)\Big) > 1$ hold. Thus, from the Lemma 3, there exist a constant $\eta > 0$, such that $r\Big(\exp(\alpha M + \beta_k\Lambda(S^*-\eta v_E)-D_k)\Big) > 1$. Furthermore, we can choose that $\eta$ and find a time $T^* > 0$ such that $|S(t)-S^*| < \eta v_E$. Hence, $\forall t_0, t > T^*$, we have,

The system $I' = (\alpha M + \beta_k\Lambda(S^*-\eta v_E) -D_k)I$ is the linear system with unique zero equilibrium. As $r\Big(\exp(\alpha M + \beta_k\Lambda(S^*-\eta v_E)-D_k)\Big) > 1$, this zero equilibrium is unstable, and the solution of that linear system will go to infinity. As we proof in Theorem 1, the sub-system constructed with $I_k$ is the cooperation system with the comparison property, which implies there exist at least one variable $I_{i, k}$ satisfies $\lim_{t\to+\infty} I_{i, k} = \infty$, which leads to the contradiction. Thus, this theorem is proved.

In the above theorem, we proved the weak persistence of the system (1.1), which is $\limsup_{t\to+\infty} I_{i, k}(t)\ge\delta^*$. But for the $\liminf_{t\to+\infty} I_{i, k}(t)\ge\delta^*$, which is also called the uniform persistence, the widely used monotone system method become noneffective. Because the relationship of different strains is competition and no cooperation between them, which is not a monotone system.

In the following theorem, we will prove the competition exclusion properties of the system (1.1).

Theorem 3 (Competition exclusion). If there exist two strains $k$ and $l$ satisfy the conditions $\beta_k = \beta_l$, and $\mu_k+\gamma_k > \mu_l+\gamma_l$, or $\beta_k < \beta_l$, and $\mu_k+\gamma_k = \mu_l+\gamma_l$, then all solutions of model (1.1) with the initial condition $(S(t_0), I_1(t_0), \cdots, I_{N_S}(t_0), R(t_0))\in A_I$, satisfy $\lim_{t\to+\infty}I_{i, k} = 0$, $i = 1, 2, \cdots, N_P$.

Proof. By solving the model (1.1), we can get

and

We first consider the conditions $\beta_k = \beta_l$, $\mu_l+\gamma_l < \mu_k+\gamma_k$, and the initial values stay in $A_I$. Thus, there exists a constant $\delta > 0$, such that $\mu_l+\gamma_l+\delta = \mu_k+\gamma_k$. Then, by $D_l+\delta E = D_k$, where $E$ is an identity matrix and the commutativity of identity matrix, $\forall I_k(t_0)$, $I_k(t)$ could be written as,

where $E$ is the identity matrix. And when $I_k(t_0) = I_l(t_0)$, the above equation could be given as,

The Lemma 2 implies there exist a constant $W > 0$, such that for all initial conditions that in $A_I$ satisfy, $0\le\|I_l(t)\|\le W$. Hence, $\lim_{t\to+\infty}\|I_k(t)\| = 0$. Thus, by the positivity of $I_k(t)$, all solutions of model (1.1) with the initial condition in $A_I$, satisfy $\lim_{t\to+\infty}I_{i, k} = 0$, $i = 1, 2, \cdots, N_P$.

When the conditions $\beta_k < \beta_l$, and $\mu_k+\gamma_k = \mu_l+\gamma_l$ hold, there exist a positive constant $\delta > 0$, such that $\beta_k+\delta = \beta_l$. Furthermore, from the boundedness of $S_i(t)$, we could find a time variable diagonal matrix $\Delta(t)$ with elements $\delta_1(t), \delta_2(t), \cdots, \delta_{N_P}(t)$, which satisfy $\forall i$, $\delta_i(t) > 0$, and $\forall i\neq j$, $0 < \delta_i(t)S_i(t) = \delta_j(t)S_j(t) = \bar{\delta} < \delta$, where $\bar{\delta} > 0$ is a constant. Thus, the right hand side of the system $\frac{d\bar{I}(t)}{dt} = (\alpha M + \beta_k\Lambda(s) + \Delta(t)\Lambda(s) -D_k)\bar{I}(t)$, satisfy $(\alpha M + \beta_k\Lambda(s) -D_k)\bar{I}(t) < (\alpha M + \beta_k\Lambda(s) + \Delta(t)\Lambda(s) -D_k)\bar{I}(t) < (\alpha M + \beta_l\Lambda(s) -D_l)\bar{I}(t)$. Then by the comparison of cooperative system and lemma 2, the solution of that system $\bar{I}(t)$ is positive and bounded.

By the commutativity of the identity matrix, then choose the same initial condition, the following equation holds.

Then similar to the above proof. Use the boundedness of $\bar{I}(t)$ and positivity of $I_k(t)$, we could prove that $\lim_{t\to+\infty}I_{i, k} = 0$, $i = 1, 2, \cdots, N_P$. Thus, we finish the proof.

Remark: The above theorem is not true for a simple assuming $R_{0, k} < R_{0, l}$. Because even for the one dimensional patch system, we could only conclude $\beta_kS^*-(\gamma_k+\mu_k+d_1) < \beta_lS^*-(\gamma_l+\mu_l+d_1)$, but not know whether the inequation $\beta_kS-(\gamma_k+\mu_k+d_1) < \beta_lS-(\gamma_l+\mu_l+d_1)$ is hold for all $S < S^*$. If this inequation doesn't hold for some $S$, the growth rate of strain $k$ may be faster than $l$, which couldn't lead to the extinction of strain $k$.

4.   Simulation

In this section, a model containing two strains and two patches is used to simulate the multiple strains' competition process in patchy environments (which is listed as model (4.1)). The baseline parameters used in the simulations of this section are $\lambda_1 = 1$, $\lambda_2 = 1.5$, $m_{21} = 0.05$, $m_{12} = 0.01$, $\mu_1 = 0.1$, $\mu_2 = 0.1$, $d_1 = 1$, $d_2 = 1$, $\gamma_1 = 1$, $\gamma_2 = 1$, $\beta_1 = 2.5$, $\beta_2 = 3$ and $\alpha = 0.8$. The baseline initial conditions are taken as $(S_1(0), S_2(0), I_{1, 1}(0), I_{1, 2}(0), I_{2, 1}(0), I_{2, 2}(0)) = (1, 2, 0.01, 0.02, 0.02, 0.0)$.

The competition exclusion property is one of the main results analyzed in the multiple strain model. In this work, the competition exclusion theorem is shown at the Theorem 3 and the simulation result that under the condition of Theorem 3 is shown in the Figure 1. The simulation results illustrate that strains with a higher infection rate could competitively exclude the other strains in a short period. This could explain why the delta strain of Covid-19 ^[20] could be the superior one in a very short period in most places of the world ^[21].

The simulation result of previous works ^[22] shows the initial condition could influence the equilibrium of the infected strain when the infection and recovery rate of different strains are the same. Therefore, in this work, we also numerical calculate the system (4.1) to verify this property. Figure 2 illustrates the simulation results of the two strains of system (4.1) are with the same infection and recovery rate. The top $2$ panels of Figure 2 illustrate the simulation result of the system (4.1) under the baseline initial condition, while the bottom $2$ panels illustrate the simulation result under the initial condition $(1, 2, 0.01, 0.02, 0.04, 0.05)$. These simulation results show that different initial conditions could change the trajectories' trend and the equilibrium of the system when the infection rate and recovery rate of that two strains were the same. These simulation results may give a biological explanation on why there are so many COVID-19 strains that could co-exist at one place and the proportion of infected cases are different in different places ^[21].

Lockdown is one strategy to slow down or stop the transmission of the epidemic. In this work, we consider the lockdown of one patch that limits or stops the immigrants from the other patches, and we call the limit the immigrants as mild lockdown and stop the immigrants is called the strict lockdown. The Figure 3 shows the simulation result of the mild (top 2 panels of Figure 3) and strict (bottom 2 panels of Figure 3) lockdown with the parameter $m_{12} = 0.001$ as the mild lockdown and the $m_{12} = 0.0$ as the strict lockdown. The other parameters are the same as the Figure 2. The simulation result implies that if the people at one patch don't change their behavior, even for a very some immigrant rate could cause the spread of the disease at this patch, and the only way to stop the transmission from one patch to another is to shut down all of the passes from this patch to the others. This means, for the global transmission of the Covid-19, unless there is no disease in all areas or everyone changes their behavior, even if all passes are mild lockdown, the spread of disease will not be prevented.

5.   Conclusions

In this work, we analyze a complex high dimensional SIR system with multiple strains and patches, and conclude the properties of that system under different basic reproductive numbers and special conditions. We find when all of the basic reproductive numbers are smaller than $1$, the disease-free equilibrium is globally stable. The persistence condition, that all basic reproductive numbers are larger than $1$, is also given. Furthermore, the competition exclusion condition for the system (1.1) is also introduced by solving the equations. By numerical analysis, we compare different lockdown strategies and find that only strict lockdown could stop the transmission of the disease in this patches.

This work gives some insights into the global dynamic of high dimension ordinary differential equations. All of the analysis methods illustrated in this work could be used in other related works about the multiple strain and patch competition. The competition exclusion theorem implies that a new strain with a higher infection rate or lower recovery rate will become the superior one, and overcome the inferior strains in the competition. This could explain why in the global spread of Covid-19, the delta strain which has a higher infection rate becomes the superior one ^[20] and the alpha strain is hard to be detected. Furthermore, according to the analysis results of this article, compared with the dominant strain delta, if the omicron strain ^[23] has stronger infection rate and similar recovery rate, then the omicron strain will replace delta strain and become the dominant one in the near future.

Acknowledgments

This work is financially supported by the GDAS' Project of Science and Technology Development (2021GDASYL-20210103089), China postdoctoral science foundation (2021M690747), National Natural Science Foundation of China (11701445, 12001139, 1201014), Science and Technology Program of Guangzhou (202007040007), GDAS' Project of Science and Technology Development (2019GDASYL-0502007), Guangdong Provincial Rural Revitalization Strategy Special Fund Project (2019KJ138), and Guangdong Basic and Applied Basic Research Foundation (2019A1515110503).