A study on the transmission dynamics of COVID-19 considering the impact of asymptomatic infection

Abstract

The COVID-19 epidemic has been spreading around the world for nearly three years, and asymptomatic infections have exacerbated the spread of the epidemic. To analyse and evaluate the role of asymptomatic infections in the spread of the epidemic, we establish an improved COVID-19 infectious disease dynamics model. We fit the epidemic data in the four time periods corresponding to the selected 614G, Alpha, Delta and Omicron variants and obtain the proportion of asymptomatic persons among the infected persons gradually increased and with the increase of the detection ratio, the cumulative number of cases has dropped significantly, but the decline in the proportion of asymptomatic infections is not obvious. Therefore, in view of the hidden transmission of asymptomatic infections, the cooperation between various epidemic prevention and control policies is required to effectively curb the spread of the epidemic.

KEYWORDS:

• COVID-19
• Asymptomatic infections
• Control reproduction number
• Lyapunov function
• Data fitting

Mathematics subjects classification 2010:

• 34A34
• 92B05

View correction statement:

Correction

1. Introduction

Since its outbreak at the end of 2019, COVID-19 has had a severe impact on countries around the world. How to effectively control the spread of the epidemic and restore normal production and life is still an issue that governments need to consider. To control the spread of COVID-19, many measures have been taken, including non-pharmacological interventions, pharmacological interventions, vaccinations, and more. However, due to the continuous emergence of variants, the epidemic prevention measures did not achieve the expected effect.

There are a large number of contagious asymptomatic infections during the spread of COVID-19 [1–4]. This feature determines the high concealment of the spread of COVID-19, and leads to an increase in the difficulty of epidemic prevention and control. There have been a large number of research results on asymptomatic infections of new coronary pneumonia, including age distribution, average proportion, and transmission intensity [5–7]. In addition, some researchers have studied the spread of COVID-19 by establishing dynamic models including asymptomatic infections. Ruan et al. established a time-varying COVID-19 transmission compartment model including asymptomatic infected persons, simulated and reviewed the development process of the Wuhan epidemic, and obtained that the asymptomatic proportion of infected persons was about 20% [8]. Rahul Subramanian et al. established a COVID-19 transmission model including asymptomatic infections, quantified asymptomatic infections in New York City, and obtained asymptomatic infections accounted for about 60% [9]. Mohamed Amouch et al. proposed a new epidemiological mathematical model of the spread of COVID-19 disease, fitted the outbreak in Monaco, and obtained that the proportion of asymptomatic patients was 30% [10]. However, these studies are all conducted for a specific variant, and cannot effectively reflect the changes caused by variant iterations.

In the more than two years since the first appearance of COVID-19, many countries and regions around the world have experienced repeated outbreaks. Taking England as an example, although multiple rounds of lockdown measures have been adopted to control the spread of the epidemic, such epidemic prevention measures have not been fully effective. When the epidemic prevention measures were gradually lifted, the epidemic rebounded again. We believe that asymptomatic infections have played a very important role in the rebound of the epidemic. The lockdown measures have effectively reduced the number of infected people to a certain extent, but because large-scale screening tests were not adopted, there were still a certain number of asymptomatic infections in the population. After the epidemic prevention measures are lifted, these undetected asymptomatic infections will cause the next round of outbreaks. If the number of asymptomatic infections can be effectively estimated in the early stage of the next outbreak, and certain epidemic prevention and control measures are taken, the spread of the epidemic can be delayed to a certain extent. We establish an improved SEIAR infectious disease dynamic model to assess the role of asymptomatic infections in the early stages of epidemic transmission.

The remainder of this article is organized as follows. In Section 2, we present the COVID-19 infectious disease model and analyse the existence and stability of equilibriums. In the third section, we use the established model to fit the actual epidemic data, and conduct a comparative analysis of the relevant kinetic parameters. In the last section, we conclude and discuss.

2. Model establishment and analysis

2.1 Model establishment

There are two typical characteristics of the spread of COVID-19. One is a large number of asymptomatic infected people with COVID-19, the other is that people with COVID-19 have a roughly 3–7 d exposure period and can show strong contagion two days before the patient shows symptoms [11,12]. Based on our knowledge of COVID-19, we have the following assumptions in our model:

(A1) We divide the exposure period into two parts: the early stage and the later stage, in which the early stage is not infectious and the later stage is infectious;

(A2) For symptomatic infections, we assume that symptomatic infections can be screened as long as they are tested, and that the detected symptomatic infections are completely isolated and no longer contagious;

(A3) Asymptomatic infections were assumed not to be tested, and deaths from illness in asymptomatic infections were not considered.

Based on the above assumptions, we establish a COVID-19 transmission compartment model including asymptomatic infected persons. The model divides the general population into susceptible ($S$), pre-exposure patients ($E_1$), late-exposure patients ($E_2$), detected symptomatic infections ($I_1$), undetected symptomatic infections ($I_2$), asymptomatic infections ($A$) and recovered ($R$). Its dynamic flow chart is shown in Figure 1.

Figure 1. Flow chart of the COVID-19 transmission dynamic model containing seven compartments. The compartments represented by the red box are infectious, and the compartments represented by the blue box are not infectious.

The corresponding propagation dynamic equation is constructed as follows: (1) $\{\begin{array}{l}{\displaystyle \frac{dS}{dt}}=\mathrm{\Lambda }-\beta S(E_2+I_2+\omega A)-\mu S,\\ {\displaystyle \frac{d{E}_1}{dt}}=\beta S(E_2+I_2+\omega A)-(\sigma +\epsilon +\mu )E_1,\\ {\displaystyle \frac{d{E}_2}{dt}}=\sigma E_1-(\alpha +\mu )E_2,\\ {\displaystyle \frac{d{I}_1}{dt}}=\rho \alpha E_2-({\gamma }_1+{\phi }_1+\mu )I_1,\\ {\displaystyle \frac{d{I}_2}{dt}}=(1-\rho )\alpha E_2-({\gamma }_2+{\phi }_2+\mu )I_2,\\ {\displaystyle \frac{dA}{dt}}=\epsilon E_1-({\gamma }_3+\mu )A,\\ {\displaystyle \frac{dR}{dt}}={\gamma }_1{I}_1+{\gamma }_2{I}_2+{\gamma }_{3}A-\mu R.\end{array}$(1) The parameters in system (1) are explained as follows. $\mathrm{\Lambda }$ is the daily number of births and $\mu$ is the daily natural mortality rate. The parameter $\beta$ represents the basal transmission rate, and $\sigma$ represents the conversion rate from patients in the pre-exposure period to patients in the late-exposure period. $\epsilon$ is the conversion rate from patients in the pre-exposure period to asymptomatic infected persons, and $\alpha$ represents the conversion rate from patients in the late exposure period to infected persons. $\omega$ is the weight of the transmission intensity of the asymptomatic infection relative to the symptomatic infection, and $\rho$ is the proportion of testing among the symptomatic infection and the value interval is [0,1]. ${\gamma }_1,{\gamma }_2,{\gamma }_3$ is the recovery rate of detected symptomatic infections, undetected symptomatic infections, and asymptomatic infections, respectively. ${\phi }_1,{\phi }_2$ is the mortality rate of detected symptomatic infections and undetected symptomatic infections, respectively. In order to maintain biological significance, the dynamic parameters above are positive real numbers.

2.2 Calculation of disease-free equilibrium and control reproduction number

The domain of the solution of the system (1) is $\mathrm{\Gamma }=\left\{(S,E_1,E_2,I_1,I_2,A,R)\in R_{+}^7|S+E_1+E_2+I_1+I_2+A+R\le \frac{\mathrm{\Lambda }}{\mu }\right\}.$Obviously, there is always a disease-free equilibrium $P_0=(S_0,0,0,0,0,0,0)$ in system (1), where $S_0=\mathrm{\Lambda }/\mu$. Then we use the next generation matrix method to calculate the control reproduction number of system (1). The Jacobian matrices $F$ and $V$ at the disease-free equilibrium are obtained from system (1) as $F=\left(\begin{array}{ccccc}0& \beta S_0& 0& \beta S_0& \omega \beta S_0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$ $V=\left(\begin{array}{ccccc}\sigma +\epsilon +\mu & 0& 0& 0& 0\\ -\sigma & \alpha +\mu & 0& 0& 0\\ 0& -\rho \alpha & {\gamma }_1+{\phi }_1+\mu & 0& 0\\ 0& -(1-\rho )\alpha & 0& {\gamma }_2+{\phi }_2+\mu & 0\\ -\epsilon & 0& 0& 0& {\gamma }_3+\mu \end{array}\right).$

The control reproduction number is the spectral radius of $F{V}^{-1}$, therefore we get. ${\mathcal{R}}_C=\frac{\beta S_0}{\sigma +\epsilon +\mu }\left(\frac{\sigma }{\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}+\frac{\epsilon \omega }{{\gamma }_3+\mu }\right).$

2.3 Existence of the endemic equilibrium

Theorem 2.3.

When ${\mathcal{R}}_C>1$, the system (1) has a unique positive equilibrium point $P^{\ast }=(S^{\ast },{E_1}^{\ast },{E_2}^{\ast },{I_1}^{\ast },{I_2}^{\ast },A^{\ast },R^{\ast })$.

Proof. Using the equilibrium equation of system (1) at the endemic equilibrium, we can get (2) $\begin{array}{rl}{S}^{\ast }& =\frac{\mathrm{\Lambda }}{\beta {E_1}^{\ast }\left(\frac{\sigma }{\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}+\frac{\epsilon \omega }{{\gamma }_3+\mu }\right)+\mu },\\ E_2^{\ast }& =\frac{\sigma }{\alpha +\mu }{E_1}^{\ast },{I_1}^{\ast }=\frac{\sigma \rho \alpha }{(\alpha +\mu )({\gamma }_1+{\phi }_1+\mu )}{E_1}^{\ast },\\ {I_2}^{\ast }& =\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}{E_1}^{\ast },A^{\ast }=\frac{\epsilon }{{\gamma }_3+\mu }{E_1}^{\ast },\\ R^{\ast }& =\frac{{E_1}^{\ast }}{\mu }\left[\frac{\sigma \rho \alpha {\gamma }_1}{(\alpha +\mu )({\gamma }_1+{\phi }_1+\mu )}+\frac{\sigma (1-\rho )\alpha {\gamma }_2}{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}+\frac{\epsilon {\gamma }_3}{{\gamma }_3+\mu }\right].\end{array}$(2)

Furthermore, according to the equilibrium equation, we obtain $\mathrm{\Lambda }-\mu S^{\ast }-(\sigma +\epsilon +\mu ){E_1}^{\ast }=0.$

Using equation (2) we can simplify the above equation to ${E_1}^{\ast }[\beta (\sigma +\epsilon +\mu ){\mathcal{R}}_C{E_1}^{\ast }+(\mu \beta S_0-\mathrm{\Lambda }\beta {\mathcal{R}}_C)]=0.$

Therefore, we can get ${E_1}^{\ast }=0,\text{or}{E_1}^{\ast }=-\frac{\mu \beta S_0-\mathrm{\Lambda }\beta {\mathcal{R}}_C}{\beta (\sigma +\epsilon +\mu ){\mathcal{R}}_C},$where $\mu \beta S_0-\mathrm{\Lambda }\beta {\mathcal{R}}_C=\mathrm{\Lambda }\beta \left(\frac{\mu S_0}{\mathrm{\Lambda }}-{\mathcal{R}}_C\right)=\mathrm{\Lambda }\beta (1-{\mathcal{R}}_C).$

Therefore, ${E_1}^{\ast }$ has a unique positive root when ${\mathcal{R}}_C>1$. According to equation (2), we can know that there is a unique positive equilibrium $P^{\ast }$ of system (1) at this time.

2.4 Stability

Theorem 2.4.1.

When ${\mathcal{R}}_C<1$, the disease-free equilibrium $P_0$ of system (1) is locally asymptotically stable.

Proof. We refer to the method of proof for stability in [13, 14]. Let the Jacobian matrix of system (1) at the disease-free equilibrium $P_0$ be $M_1$. The matrix $M_1$ is detailed in the Appendix A. When ${\mathcal{R}}_C<1$, assume that there is ${\lambda }_0$ satisfied $\mathrm{Re}({\lambda }_0)\ge 0$, and $det({\lambda }_{0}E-M_1)=0$. Expanding the characteristic equation, we can obtain. (3) $\begin{array}{rl}& ({\lambda }_0+\mu {)}^2({\lambda }_0+\alpha +\mu )({\lambda }_0+{\gamma }_1+{\phi }_1+\mu )({\lambda }_0+{\gamma }_2+{\phi }_2+\mu )({\lambda }_0+{\gamma }_3+\mu )\\ & \cdot \{{\lambda }_0+(\sigma +\epsilon +\mu )[1-F({\lambda }_0)]\}=0,\end{array}$(3) Where $\begin{array}{rl}F({\lambda }_0)& =\frac{\beta S_0}{\sigma +\epsilon +\mu }[\frac{\sigma }{{\lambda }_0+\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{({\lambda }_0+\alpha +\mu )({\lambda }_0+{\gamma }_2+{\phi }_2+\mu )}\\ & +\frac{\epsilon \omega }{{\lambda }_0+{\gamma }_3+\mu }].\end{array}$

Since $\mathrm{Re}({\lambda }_0)\ge 0$ is assumed, the characteristic equation (3) can be simplified to. ${\lambda }_0+(\sigma +\epsilon +\mu )[1-F({\lambda }_0)]=0.$

The necessary condition for the above formula to be established is (4) $\mathrm{Re}[F({\lambda }_0)]\ge 1.$(4)

Furthermore, assumed $\mathrm{Re}({\lambda }_0)\ge 0$, we can obtain. (5) $\begin{array}{rl}|F({\lambda }_0)|& =|\frac{\beta S_0}{\sigma +\epsilon +\mu }[\frac{\sigma }{{\lambda }_0+\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{({\lambda }_0+\alpha +\mu )({\lambda }_0+{\gamma }_2+{\phi }_2+\mu )}\\ & +\frac{\epsilon \omega }{{\lambda }_0+{\gamma }_3+\mu }]|\\ & \le \frac{\beta S_0}{\sigma +\epsilon +\mu }[\left|\frac{\sigma }{{\lambda }_0+\alpha +\mu }\right|+\left|\frac{\sigma (1-\rho )\alpha }{({\lambda }_0+\alpha +\mu )({\lambda }_0+{\gamma }_2+{\phi }_2+\mu )}\right|\\ & +\left|\frac{\epsilon \omega }{{\lambda }_0+{\gamma }_3+\mu }\right|]\\ & \le \frac{\beta S_0}{\sigma +\epsilon +\mu }\left[\left|\frac{\sigma }{\alpha +\mu }\right|+\left|\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}\right|+\left|\frac{\epsilon \omega }{{\gamma }_3+\mu }\right|\right]\\ & =\frac{\beta S_0}{\sigma +\epsilon +\mu }\left[\frac{\sigma }{\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}+\frac{\epsilon \omega }{{\gamma }_3+\mu }\right]\\ & ={\mathcal{R}}_C.\end{array}$(5)

Therefore, in the case of ${\mathcal{R}}_C<1$, it can be known from inequality (5) that (6) $\mathrm{Re}[F({\lambda }_0)]<1.$(6)

At this point, it can be seen that there is a contradiction between inequalities (4) and (6), so the initial assumption about ${\lambda }_0$ is invalid, that is, all the characteristic roots of the characteristic equation $det({\lambda }_{0}E-M_1)=0$ have negative real parts. That is to say, when ${\mathcal{R}}_C<1$, the disease-free equilibrium $P_0$ of system (1) is locally asymptotically stable.

Theorem 2.4.2.

When ${\mathcal{R}}_C<1$, the disease-free equilibrium $P_0$ of system (1) is globally asymptotically stable.

Proof. Define function $h(x)=x-1-\ln x,x>0.$For any $x>0$ there is $h(x)\ge 0$, and $h(x)=0$ if and only if $x=1$.

Let $U(t)=(S(t),E_1(t),E_2(t),I_1(t),I_2(t),A(t),R(t))$ be the solution of system (1) through any $\varphi :=({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4,{\varphi }_5,{\varphi }_6,{\varphi }_7)\in {\mathbb{R}}_{+}^7$. Clearly, $U(t)$ is bounded and $S(t)>0$ holds for any $t\ge 0$. Referring to the construction method in [15], define a function $V$ on $\mathrm{\Omega }=\{\varphi \in R_{+}^7:{\varphi }_1>0\}$ as follows, $\begin{array}{l}V(\varphi )=S_{0}h\left(\frac{{\varphi }_1}{S_0}\right)+{\mathcal{R}}_C{\varphi }_2+\frac{\beta S_0}{\alpha +\mu }{\varphi }_3+\frac{\omega \beta S_0}{{\gamma }_3+\mu }{\varphi }_6\\ +\frac{\beta S_0}{{\gamma }_2+{\phi }_2+\mu }[(1-\rho ){\varphi }_3+{\varphi }_5]-\frac{\beta S_0(1-\rho )\mu }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}{\varphi }_3.\end{array}$

The derivative of $V$ with respect to $t(t>0)$ along $U(t)$ is $\begin{array}{rl}\frac{dV(U(t))}{dt}& =\left(1-\frac{S_0}{S(t)}\right)[\mathrm{\Lambda }-\beta S(t)(E_2(t)+I_2(t)+\omega A(t))-\mu S(t)]\\ & +{\mathcal{R}}_C[\beta S(t)(E_2(t)+I_2(t)+\omega A(t))-(\sigma +\epsilon +\mu )E_1(t)]\\ & +\frac{\beta S_0}{\alpha +\mu }[\sigma E_1(t)-(\alpha +\mu )E_2(t)]+\frac{\omega \beta S_0}{{\gamma }_3+\mu }[\epsilon E_1(t)-({\gamma }_3+\mu )A(t)]\\ & +\frac{\beta S_0}{{\gamma }_2+{\phi }_2+\mu }[\sigma (1-\rho )E_1(t)-\mu (1-\rho )E_2(t)-({\gamma }_2+{\phi }_2+\mu )I_2(t)]\\ & -\frac{\beta S_0(1-\rho )\mu }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}[\sigma E_1(t)-(\alpha +\mu )E_2(t)]\\ & =-\frac{\mu }{S(t)}[S(t)-S_0{]}^2+({\mathcal{R}}_C-1)\beta S(t)(E_2(t)+I_2(t)+\omega A(t)).\end{array}$

Hence when ${\mathcal{R}}_C<1$, we have $(dV(t)/dt)\le 0$, which means that $V$ is a Lyapunov function on $\{U(t):t\ge 1\}\in \mathrm{\Omega }$. And $dV(t)/dt=0$ if and only if $U(t)=P_0$. Therefore, the largest invariant set of $\{U(t)|dV(t)/dt=0,U(t)\in \mathrm{\Omega }\}$ is $\{P_0\}$. According to the LaSalle's Invariance Principle, when ${\mathcal{R}}_C<1$, the disease-free equilibrium $P_0$ of system (1) is globally asymptotically stable.

Theorem 2.4.3.

When ${\mathcal{R}}_C>1$, the disease-free equilibrium $P_0$ of system (1) is unstable.

Proof. Similar to the proof method in Theorem 2.4.1, when ${\mathcal{R}}_C>1$, assumed that there is ${\lambda }_1$ satisfied $\mathrm{Re}({\lambda }_1)\ge 0$, and $det({\lambda }_{1}E-M_1)=0$. Expanding the characteristic equation, we can obtain. (7) $\begin{array}{rl}& ({\lambda }_1+\mu {)}^2({\lambda }_1+\alpha +\mu )({\lambda }_1+{\gamma }_1+{\phi }_1+\mu )({\lambda }_1+{\gamma }_2+{\phi }_2+\mu )({\lambda }_1+{\gamma }_3+\mu )\\ & \cdot \{{\lambda }_1+(\sigma +\epsilon +\mu )[1-F({\lambda }_1)]\}=0,\end{array}$(7) where $\begin{array}{rl}F({\lambda }_1)& =\frac{\beta S_0}{\sigma +\epsilon +\mu }[\frac{\sigma }{{\lambda }_1+\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{({\lambda }_1+\alpha +\mu )({\lambda }_1+{\gamma }_2+{\phi }_2+\mu )}\\ & +\frac{\epsilon \omega }{{\lambda }_1+{\gamma }_3+\mu }].\end{array}$

The characteristic equation (7) can be simplified to (8) $({\lambda }_1+\alpha +\mu )({\lambda }_1+{\gamma }_2+{\phi }_2+\mu )({\lambda }_1+{\gamma }_3+\mu )\{{\lambda }_1+(\sigma +\epsilon +\mu )[1-F({\lambda }_1)]\}=0.$(8)

To simplify the subsequent formulation, redefine the notation as $\begin{array}{l}{B}_1=\alpha +\mu ,B_2={\gamma }_2+{\phi }_2+\mu ,B_3={\gamma }_3+\mu ,C_1=\sigma +\epsilon +\mu ,\\ C_2=\sigma ,C_3=\sigma (1-\rho )\alpha ,C_4=\epsilon \omega ,D=\beta S_0.\end{array}$

Equation (8) is now rewritten as. $\begin{array}{rl}& ({\lambda }_1+B_1)({\lambda }_1+B_2)({\lambda }_1+B_3)\\ & \times \left\{{\lambda }_1+C_1\left\{1-\frac{D}{C_1}\left[\frac{C_2}{{\lambda }_1+B_1}+\frac{C_3}{({\lambda }_1+B_1)({\lambda }_1+B_2)}+\frac{C_4}{{\lambda }_1+B_3}\right]\right\}\right\}=0.\end{array}$

The result of expanding the above formula is ${{\lambda }_1}^4+a_1{{\lambda }_1}^3+a_2{{\lambda }_1}^2+a_3{\lambda }_1+a_4=0.$

Where $\begin{array}{rl}{a}_1& =B_1+B_2+B_3+C_1,\\ a_2& =B_1{B}_2+B_1{B}_3+B_2{B}_3+B_1{C}_1+B_2{C}_1+B_3{C}_1-D{C}_2-D{C}_4,\\ a_3& =B_1{B}_2{B}_3+B_1{B}_2{C}_1+B_1{B}_3{C}_1+B_2{B}_3{C}_1-D{B}_2{C}_2-D{B}_3{C}_2\\ & -D{C}_3-D{B}_1{C}_4-D{B}_2{C}_4,\\ a_4& =B_1{B}_2{B}_3{C}_1-D{B}_2{B}_3{C}_2-D{B}_1{B}_2{C}_4-D{B}_3{C}_3.\end{array}$

Define function $G(\lambda )={\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4.$

Under the new notation, the control reproduction number is ${\mathcal{R}}_C=\frac{D}{C_1}\left(\frac{C_2}{B_1}+\frac{C_3}{B_1{B}_2}+\frac{C_4}{B_3}\right).$Therefore, when ${\mathcal{R}}_C>1$, there is $G(0)=B_1{B}_2{B}_3{C}_1-D{B}_2{B}_3{C}_2-D{B}_1{B}_2{C}_4-D{B}_3{C}_3=\frac{1}{B_1{B}_2{B}_3{C}_1}(1-{\mathcal{R}}_C)<0.$

In addition, there is $G(+\mathrm{\infty })=+\mathrm{\infty }$, so equation (8) has at least one positive root. That is to say, the characteristic equation $det({\lambda }_{1}E-M_1)=0$ has at least one positive root, that is, the assumption is true. This indicates that, when ${\mathcal{R}}_C>1$, the disease-free equilibrium $P_0$ of system (1) is unstable.

Theorem 2.4.4.

When ${\mathcal{R}}_C>1$, the system (1) is locally asymptotically stable at the endemic equilibrium $P^{\ast }$.

Proof. Let the Jacobian matrix of system (1) at the endemic equilibrium $P^{\ast }$ be $M_2$. The matrix $M_2$ is detailed in the Appendix B. When ${\mathcal{R}}_C>1$, assume that there is ${\lambda }_2$ satisfying $\mathrm{Re}({\lambda }_2)\ge 0$, such that $det({\lambda }_{2}E-M_2)=0$. Expanding the characteristic equation, we can obtain. (9) $\begin{array}{rl}& ({\lambda }_2+\mu )({\lambda }_2+\alpha +\mu )({\lambda }_2+{\gamma }_1+{\phi }_1+\mu )({\lambda }_2+{\gamma }_2+{\phi }_2+\mu )({\lambda }_2+{\gamma }_3+\mu )\\ & \cdot \{({\lambda }_2+\mu )\left[{\lambda }_2+(\sigma +\epsilon +\mu )\left(1-\frac{S^{\ast }}{S_0}F({\lambda }_2)\right)\right]\\ & +\beta ({E_2}^{\ast }+{I_2}^{\ast }+\omega A^{\ast })({\lambda }_2+\sigma +\epsilon +\mu )\phantom{\left[{\lambda }_2+(\sigma +\epsilon +\mu )\left(1-\frac{S^{\ast }}{S_0}F({\lambda }_2)\right)\right]}\}=0,\end{array}$(9) where, $\begin{array}{rl}F({\lambda }_2)& =\frac{\beta S_0}{\sigma +\epsilon +\mu }[\frac{\sigma }{{\lambda }_2+\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{({\lambda }_2+\alpha +\mu )({\lambda }_2+{\gamma }_2+{\phi }_2+\mu )}\\ & +\frac{\epsilon \omega }{{\lambda }_2+{\gamma }_3+\mu }].\end{array}$

The characteristic equation (9) can be simplified to $\begin{array}{rl}& ({\lambda }_2+\mu )({\lambda }_2+\sigma +\epsilon +\mu )+\beta ({E_2}^{\ast }+{I_2}^{\ast }+\omega A^{\ast })({\lambda }_2+\sigma +\epsilon +\mu )\\ & =\frac{S^{\ast }(\sigma +\epsilon +\mu )}{S_0}({\lambda }_2+\mu )F({\lambda }_2).\end{array}$

Using equation (2), the above equation can be transformed into $\begin{array}{rl}& ({\lambda }_2+\mu )({\lambda }_2+\sigma +\epsilon +\mu )+\frac{{\mathcal{R}}_C{E_1}^{\ast }}{S_0}(\sigma +\epsilon +\mu )({\lambda }_2+\sigma +\epsilon +\mu )\\ & =\frac{S^{\ast }(\sigma +\epsilon +\mu )}{S_0}({\lambda }_2+\mu )F({\lambda }_2).\end{array}$

The above equation can be simplified to (10) $1+\frac{{\mathcal{R}}_C{E_1}^{\ast }(\sigma +\epsilon +\mu )}{S_0}\cdot \frac{1}{{\lambda }_2+\mu }=\frac{S^{\ast }(\sigma +\epsilon +\mu )}{S_0({\lambda }_2+\sigma +\epsilon +\mu )}F({\lambda }_2).$(10)

Furthermore, from the equilibrium equation we can obtain $\beta S^{\ast }({E_2}^{\ast }+{I_2}^{\ast }+\omega A^{\ast })=(\sigma +\epsilon +\mu ){E_1}^{\ast }.$

Using equation (2) to replace the parameters of the above equation, we can obtain $\beta S^{\ast }\left(\frac{\sigma }{\alpha +\mu }+\frac{\sigma (1-\rho )\alpha }{(\alpha +\mu )({\gamma }_2+{\phi }_2+\mu )}+\frac{\omega \epsilon }{{\gamma }_3+\mu }\right){E_1}^{\ast }=(\sigma +\epsilon +\mu ){E_1}^{\ast },$and ${\mathcal{R}}_C=\frac{S_0}{S^{\ast }}.$

Similar to the proof process of inequality (5), we get $|F({\lambda }_2)|\le {\mathcal{R}}_C,$taking the norm on the right-hand side of equation (10), we get (11) $\begin{array}{r}\left|\frac{S^{\ast }(\sigma +\epsilon +\mu )}{S_0({\lambda }_2+\sigma +\epsilon +\mu )}F({\lambda }_2)\right|\le \left|\frac{S^{\ast }(\sigma +\epsilon +\mu )}{S_0({\lambda }_2+\sigma +\epsilon +\mu )}{\mathcal{R}}_C\right|=\left|\frac{\sigma +\epsilon +\mu }{{\lambda }_2+\sigma +\epsilon +\mu }\right|\le 1.\end{array}$(11)

Since we assume $\mathrm{Re}({\lambda }_2)\ge 0$, there is $\mathrm{Re}[1/({\lambda }_2+\mu )]>0$, and the left-hand side of equation (10) satisfies $\mathrm{Re}\left(1+\frac{{\mathcal{R}}_C{E_1}^{\ast }(\sigma +\epsilon +\mu )}{S_0}\cdot \frac{1}{{\lambda }_2+\mu }\right)=1+\frac{{\mathcal{R}}_C{E_1}^{\ast }(\sigma +\epsilon +\mu )}{S_0}\mathrm{Re}\left(\frac{1}{{\lambda }_2+\mu }\right)>1.$

Therefore, taking the norm on the left-hand side of equation (10), we get (12) $\left|1+\frac{{\mathcal{R}}_C{E_1}^{\ast }(\sigma +\epsilon +\mu )}{S_0}\cdot \frac{1}{{\lambda }_2+\mu }\right|>1.$(12)

At this point we can see that the inequalities (11) and (12) are contradictory to the equation (10), so the initial assumption about ${\lambda }_2$ is invalid, that is, all the eigenvalues of the characteristic equation $det({\lambda }_{2}E-M_2)=0$ have negative real parts. That is, when ${\mathcal{R}}_C>1$, the system (1) is locally asymptotically stable at the endemic equilibrium $P^{\ast }$.

3. Numerical fitting

3.1 Selection of fitting objects

The SARS-CoV-2 has mutated several times since it was first discovered in late 2019. Among them, there are 7 variants of concern that caused a global pandemic, including 614D, 614G, Alpha, Beta, Gamma, Delta and Omicron variants [16]. 614D is an original variant and there is a lack of data on its infectors; although the Beta variant has strong immune evasion ability and its transmission is stronger than that of the Alpha variant, its epidemic process is quickly overtaken by the more transmissible variant; similarly, the Gamma variant, although more pathogenic, was rapidly replaced by the Delta variant in its epidemic process due to its weaker transmission; therefore, we do not consider the above three variants. This study focuses only on four of the major virus variants that are widely circulating worldwide, namely 614G, Alpha, Delta and Omicron variants.

For the 614G variant, we select the epidemic data from 1 June 2020 to 31 July 2020 in England for simulation, which was recorded as the first time period. Since then, the spread of COVID-19 in England has slowed until early November 2020, when a second nationwide lockdown began and lasted for about a month. The second national lockdown did play a relatively obvious role in curbing the spread of the epidemic, but due to the dual effects of the gradual relaxation of the lockdown policy and the emergence of the Alpha variant, England has experienced the third large-scale outbreak in a short period of time [17]. Therefore, for the Alpha variant, we select the epidemic data from 1 December 2020 to 31 December 2020 in England for simulation, which was recorded as the second time period.

Since then, England has entered the third national lockdown and lasted for about two months, which has brought the epidemic in England under certain control. However, similarly, England began its third large-scale outbreak in May 2021 due to the dual effects of the gradual easing of the third lockdown policy and the emergence of the Delta variant. Since the Delta variant first appeared in England in mid-April 2021, it has replaced the Alpha variant as the predominant variant in England in just over a month. Therefore, for the Delta variant, we select the epidemic data from 1 June 2021 to 31 June 2021 in England for simulation, which is recorded as the third time period. Since then, the UK has been adopting a lax epidemic prevention policy, which has kept the daily new case data at a high position.

For the Omicron variant, we select the epidemic data from 18 March 2022 to 9 April 2022 in Shanghai, China for simulation, which is recorded as the fourth time period [18]. There are two main reasons for choosing the epidemic data in Shanghai during this period as the research object. First, this round of epidemic in Shanghai is basically caused by single Omicron variant, which facilitates our comparative analysis. Secondly, this time period is the early stage of the outbreak in Shanghai, which is more suitable for the mathematical model. The relative time stage of the fitting time period selected relative to this round of epidemics are shown in Figure 2. The first three fitting time periods we selected are the early stages of three large-scale outbreaks in England, and the fourth time period is the early stages of the Shanghai outbreak in April 2022. We use the established model to fit the epidemic data in these four time periods, and comparatively analyse the role of asymptomatic infections in the spread of the epidemic.

Figure 2. Relative position map of the four data fitting time periods during the corresponding outbreak. a, b and c, The relative position map of the selected first, second and third data fitting time periods during the first three outbreaks of the corresponding epidemic in England. d, The relative position map of the selected fourth data fitting time period during the corresponding outbreak in Shanghai in April 2022.

3.2 Numerical fitting results

For the selected 614G, Alpha, Delta and Omicron variants, based on the established model, we use the least squares method to simulate the epidemic data in different periods. Among them, the epidemic data corresponding to the 614G, Alpha and Delta variants are all from the official website of the British government (https://coronavirus.data.gov.uk/), and the epidemic data corresponding to the Omicron variant are from the Shanghai Municipal Health Commission (https://wsjkw.sh.gov.cn/).

The fitting results of the epidemic data in the four time periods are shown in Figure 3, in which Figures 3(a–d) show the data fitting in the first to fourth time periods in turn. From the figure, we can see that the simulation results of the epidemic data in the four time periods are relatively satisfactory, which indicates that our model has good universality. Due to the small value and large fluctuation of the original data, the simulation results of the first time period and the fourth time period are not perfect in some parts, but they are generally acceptable. The related dynamic parameters obtained by fitting are shown in Table 1.

Figure 3. Simulation results of the daily number of newly confirmed cases of epidemic data in four time periods. a, Simulation results in the first time period, corresponding to the 614G variant. b, Simulation results in the second time period, corresponding to the Alpha variant. c, Simulation results in the third time period, corresponding to the Delta variant. d, Simulation results in the fourth time period, corresponding to the Omicron variant.

Table 1. Fitting results of dynamic parameters in model (1).

Parameter	Value 614G	Alpha	Delta	Omicron	Reference
Λ	1740	1740	1740	216	[19,20]
μ	2.5753*10^∧(−5)	2.5753*10^∧(−5)	2.5753*10^∧(−5)	2.4303*10^∧(−5)	[19,20]
β	5.3720*10^∧(−9)	7.2151*10^∧(−9)	9.0205*10^∧(−9)	3.2493*10^∧(−8)	Fitted
σ	0.1975	0.1975	0.1975	0.1975	[21]
ε	0.3415	0.5748	0.6768	0.5745	Fitted
α	0.5000	0.5000	0.5000	0.5000	[22]
ω	0.6524	0.6524	0.6524	0.6524	[7]
ρ	0.4689	0.1103	0.2266	0.5266	Fitted
γ1	0.0588	0.0588	0.0588	0.0588	[23,24,25]
γ2	0.0769	0.0811	0.0704	0.0537	Fitted
γ3	0.2770	0.3746	0.4810	0.4149	Fitted
φ1	1.782^6*10^(−5)	1.7826*10^(−5)	1.7826*10^(−5)	1.7826*10^(−5)	[26]
φ2	5.5963*10^(−3)	4.4054*10^(−3)	5.0410*10^(−3)	5.0179*10^(−3)	Fitted

In addition, we performed a sensitivity analysis of the system dynamic parameters. Sensitivity analysis can be used to obtain the effect of parameters on the system and help to evaluate control measures. The sensitivity of the control reproduction number ${\mathcal{R}}_C$ with respect to any parameter $x$ is calculated as ${\xi }_{{\mathcal{R}}_C}^x=\frac{\mathrm{\partial }{\mathcal{R}}_C}{\mathrm{\partial }x}\cdot \frac{x}{{\mathcal{R}}_C}.$

Using the first set of parameter values obtained by data fitting, the sensitivity of the control reproduction number with respect to the corresponding parameters is obtained, as shown in Figure 4. From the figure, we can see that there are four parameters with elasticity values greater than 0, which indicates that the spread of the epidemic will expand when these parameters increase, and the spread of the epidemic will shrink when the remaining parameters with elasticity values less than 0 increase. We can also see the impact of control measures on the epidemic. Taking parameters with elastic values less than 0 as an example, such as the proportion of testing among the symptomatic infection ($\rho$), the recovery rate of undetected symptomatic infections (${\gamma }_2$), and the recovery rate of asymptomatic infections (${\gamma }_3$), as long as the values of these three parameters are increased by control measures, the spread of the epidemic can be effectively delayed.

Figure 4. The sensitivity of the control reproduction number with respect to the corresponding parameters.

3.3 Horizontal comparison of fitting results

From the perspective of mathematical model simulation, we verify the characteristics of enhanced infectivity and weakened symptoms of infected persons of the four SARS-CoV-2 variants concerned in this study, as shown in Figure 5. We count the proportions of the fitting results of the three types of infected persons. From the circular distribution map of the types of infected people, we can see that the proportion of asymptomatic patients in the four time periods was 38.05%, 52.89%, 54.02% and 62.03% respectively, showing a clear upward trend. This result shows that the overall symptom level of the infected person is decreasing, which is consistent with the clinical statistics of the decrease in the severe rate of the infected person.

Figure 5. Horizontal comparison results of epidemic data simulation in four time periods. a, b, d and e, Changes in the proportion of the three types of infected persons. c, Model prediction results of the peak time of the epidemic in four time periods. f, Model predictions of maximum daily confirmed cases over four time periods.

In addition, we also compare the propagation speed and propagation strength of different variants. From the perspective of the model, we predict the subsequent development of the epidemic in four time periods. From Figures 5(c and f), we can see that the peak time of the epidemic development process caused by the variants of interest in this study gradually shortens, and the peak value of daily confirmed cases gradually increases, and there is even an order of magnitude difference. This indicates that the speed of transmission, as well as the intensity of transmission, of the four COVID-19 variants of interest in this study has increased. In addition, it can also be seen from the dynamic parameters obtained by fitting that the basic propagation rate has increased significantly.

3.4 The impact of the rate of participation in testing among symptomatic infected persons on the spread of the epidemic

Although the rate of testing among those with symptomatic infection is fixed in the model simulation, in practice, with the increase in the number of newly diagnosed patients every day, certain prevention and control measures will inevitably be taken to curb the spread of the epidemic. The most direct and effective prevention and control measure is to increase the detection intensity, which is reflected in the model as the increase of the parameter $\rho$. Therefore, we explore the effect of the change of parameter $\rho$ on the spread of the epidemic.

Based on the simulation results over four time periods, while keeping other parameters unchanged, re-prediction is carried out for the cases of $\rho$ = 20%, 40%, 60%, and 80%, respectively, as shown in Figures 6 and 7. From the results, we can see that with the increase of $\rho$, the cumulative number of infected people has dropped significantly. In addition, the proportion of asymptomatic infections has declined to a certain extent, but the decline is not significant. As the proportion of people tested among symptomatic infections increases, the impact of symptomatic infection on the spread of the epidemic gradually decreases. Since the model assumes that asymptomatic infected people will not participate in the detection, although the transmission intensity is lower than that of symptomatic infected people, it will still cause a certain amount of transmission. Therefore, even if the proportion of symptomatic infected people participating in the test reaches 80%, there are still a certain proportion of asymptomatic infections.

Figure 6. The impact of the tested ratio ρ among symptomatic infected persons on the spread of the epidemic. a, The influence of parameter ρ changes on the cumulative amount of infected people, based on the fitting results of the first time period. b, The impact of parameter ρ changes on the proportion of different types of infected people, based on the fitting results of the first time period. c, The influence of parameter ρ changes on the cumulative amount of infected people, based on the fitting results of the second time period. d, The impact of parameter ρ changes on the proportion of different types of infected people, based on the fitting results of the second time period.

Figure 7. The impact of the tested ratio ρ among symptomatic infected persons on the spread of the epidemic. a, The influence of parameter ρ changes on the cumulative amount of infected people, based on the fitting results of the third time period. b, The impact of parameter ρ changes on the proportion of different types of infected people, based on the fitting results of the third time period. c, The influence of parameter ρ changes on the cumulative amount of infected people, based on the fitting results of the fourth time period. d, The impact of parameter ρ changes on the proportion of different types of infected people, based on the fitting results of the fourth time period.

In addition, we compare the differences between the different variants when the rate of detection in symptomatic infected individuals changed. As can be seen from Figure 8, when the detection ratio of symptomatic infections increases from 20% to 80%, the decline ratio of the total number of infected people decreases slowly with the variation of variants. However, the decline rates remain above 70%, indicating that increased testing efforts can still effectively reduce the total number of infected people. In addition, the proportion of asymptomatic infections decrease significantly with the change of variants. From Figure 8, we can see that for the Omicron variant, when the detection ratio increased from 20% to 80%, the proportion of asymptomatic infections drop by only about 10%. The specific data are shown in Table 2, where ${\zeta }_1$ represents the percentage of decline in the total number of infected people, and ${\zeta }_2$ represents the percentage of decline in the number of asymptomatic infections. From the results of this part, we can see that in the process of epidemic prevention, only improving the detection of symptomatic infections cannot completely block the spread of the epidemic. It is necessary to cooperate with certain social isolation policies to effectively block hidden transmission caused by asymptomatic infections.

Figure 8. The impact of the tested ratio ρ among symptomatic infected persons on the spread of the epidemic. a, The reduction rate in the total number of infected people when the detection rate goes from 20% to 80%. b, The reduction rate in asymptomatic infections when testing rates go from 20% to 80%.

Table 2. The impact of parameter ρ changes on the spread of the epidemic in different time periods.

	614G	Alpha	Delta	Omicron
${\zeta }_1$	96.51%	78.37%	83.31%	77.30%
${\zeta }_2$	40.90%	35.37%	34.66%	11.72%

3.5 Model-based prediction results

We use the established model to predict the spread of the epidemic. Since the English government no longer publishes epidemic data from May 2022, we choose WHO Coronavirus (COVID-19) Dashboard (https://covid19.who.int/) as the source of epidemic data to fit and predict the development of the epidemic in the UK. As shown in Figure 9, we forecast the development of the epidemic in the UK from late 2022 to early 2023. Predictions suggest a larger outbreak in the UK in the winter of 2022, which would last for around four months. In order to avoid the collapse of the healthcare system and the negative impact on productive life, the UK should take the necessary measures to hinder the spread of the epidemic.

Figure 9. The latest UK epidemic fitting and forecast results. Based on the predicted results, the outbreak will last from September 2022 to January 2023. The number of confirmed cases will peak in early November 2022.

4. Conclusion

In this study, a SEIAR infectious disease dynamic model was established. We first analyse the local stability at the equilibrium point of the model system. As a result, when ${\mathcal{R}}_C<1$, the disease-free equilibrium $P_0$ of the system is locally asymptotically stable; when ${\mathcal{R}}_C>1$, the disease-free equilibrium $P_0$ of the system is unstable and the endemic equilibrium $P^{\ast }$ is locally asymptotically stable. In addition, we use the method of constructing the Lyapunov function to prove that the disease-free equilibrium $P_0$ of the system is globally asymptotically stable when ${\mathcal{R}}_C<1$. In addition, we use the model of the architecture to numerically fit the COVID-19 epidemic data. We select corresponding data for the four main variants, 614G, Alpha, Delta and Omicron, and performed fitting analysis. The comparative analysis of the fitting results shows that the proportion of asymptomatic infections among infected persons is gradually increasing, and the propagation speed and propagation strength of the variant are increased to a certain extent. Finally, we analyse the impact of the detection ratio of symptomatic infections on the spread of the epidemic. With the increase of the detection ratio, the cumulative cases drop significantly, but asymptomatic infections still maintain a certain proportion.

5. Discussion

In this paper, the epidemic data of the main epidemic variants of the novel coronavirus are integrated, and the proportion of asymptomatic infected people during the epidemic of different variants is analysed by mathematical model method, and a horizontal comparative analysis is carried out [27].

Since asymptomatic infections are difficult to be actively detected, the transmission caused by asymptomatic infections is difficult to be interrupted under relatively loose epidemic prevention policies. Especially after the end of the last round of the epidemic, the epidemic prevention policy tends to be gradually loosened as the number of daily confirmed cases declines, which will cause the number of asymptomatic infections to accumulate again. It can also be seen from the comparison results of this study that in general, the pathogenicity of the SARS-CoV-2 variants is decreasing but the infectivity is increasing, which means that the role of asymptomatic infected persons in the spread of the epidemic has increased significantly. And due to increased infectivity, the absolute number of symptomatic infections will not necessarily decrease. This shows that epidemic prevention and control measures are still very important, otherwise there will be continuous challenges to the medical carrying capacity.

Large-scale social activities will also cause large-scale spread of the epidemic in the short term. Taking the Tokyo Olympics as an example, during the entire event period, not only the cities where the competition were held, but also other cities that did not hold the competition experienced outbreaks [28]. During the Tokyo Olympics, Japan's epidemic prevention policy was relaxed, and social mobility increased. This made a large number of asymptomatic infections accumulated in the crowd, resulting in a full-scale outbreak of the epidemic. Therefore, when holding large-scale social activities, it is more necessary to strengthen the control of social mobility to reduce the probability of large-scale outbreaks in the short term.

In view of the increasing role of asymptomatic infections in the spread of the epidemic, non-pharmaceutical intervention (NPI) may play a more important role [29]. Common NPI measures include wearing masks, maintaining social distance, etc. Compared with symptomatic infected individuals, asymptomatic infected individuals have lower viral load and weaker transmission ability, which enables NPI to better block the transmission caused by asymptomatic infected individuals. In addition, vaccination policies still need to be actively implemented, as vaccination is the most effective way to reduce the risk of transmission overall. Therefore, in the process of epidemic prevention and control, when the detection capacity is limited, NPI and vaccination are important measures and means to avoid the collapse of the medical system.

This paper also has some limitations. The model in this paper is only suitable for numerical fitting at the beginning of the outbreak, and does not consider more factors affecting the epidemic. In addition, reinfection due to decreased antibody levels in recovered people is not included in the model, so the model is only suitable for short-term outbreak data fitting.

Acknowledgments

We thank all the individuals who generously shared their time and materials for this study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article was originally published with errors, which have now been corrected in the online version. Please see Correction (http://dx.doi.org/10.1080/17513758.2023.2259761)

Additional information

Funding

This study was funded by the Development and Reform Commission of Guangdong Province [grant number 2017KTSCX124] and Natural Science Foundation of China (NSFC 11901027); Postgraduate Teaching Research and Quality Improvement Project of BUCEA; Graduate Research and Innovation Project of BUCEA.

Appendices

Appendix A

Jacobian matrix of system (1) $M_1$ at the disease-free equilibrium $P_0$. $\begin{array}{rl}{M}_1& =(\begin{array}{cccc}-\mu & 0& -\beta S_0& 0\\ 0& -(\sigma +\epsilon +\mu )& \beta S_0& 0\\ 0& \sigma & -(\alpha +\mu )& 0\\ 0& 0& \rho \alpha & -({\gamma }_1+{\phi }_1+\mu )\\ 0& 0& (1-\rho )\alpha & 0\\ 0& \epsilon & 0& 0\\ 0& 0& 0& {\gamma }_1\end{array}\\ & \begin{array}{ccc}\beta S_0& -\omega \beta S_0& 0\\ \beta S_0& \omega \beta S_0& 0\\ 0& 0& 0\\ 0& 0& 0\\ -({\gamma }_2+{\phi }_2+\mu )& 0& 0\\ 0& -({\gamma }_3+\mu )& 0\\ {\gamma }_2& {\gamma }_3& -\mu \end{array})\end{array}$

Appendix B

The Jacobian matrix of system (1) $M_2$ at the endemic equilibrium $P^{\ast }$. $\begin{array}{rl}{M}_2& =(\begin{array}{cccc}-\beta ({E_2}^{\ast }+{I_2}^{\ast }+\omega A^{\ast })-\mu & 0& -\beta S^{\ast }& 0\\ \beta ({E_2}^{\ast }+{I_2}^{\ast }+\omega A^{\ast })& -(\sigma +\epsilon +\mu )& \beta S^{\ast }& 0\\ 0& \sigma & -(\alpha +\mu )& 0\\ 0& 0& \rho \alpha & -({\gamma }_1+{\phi }_1+\mu )\\ 0& 0& (1-\rho )\alpha & 0\\ 0& \epsilon & 0& 0\\ 0& 0& 0& {\gamma }_1\end{array}\\ & \begin{array}{ccc}-\beta S^{\ast }& -\omega \beta S^{\ast }& 0\\ 0& \beta S^{\ast }& \omega \beta S^{\ast }\\ 0& 0& 0\\ 0& 0& 0\\ -({\gamma }_2+{\phi }_2+\mu )& 0& 0\\ 0& -({\gamma }_3+\mu )& 0\\ {\gamma }_2& {\gamma }_3& -\mu \end{array})\end{array}$