A COVID-19 Infection Model Considering the Factors of Environmental Vectors and Re-Positives and Its Application to Data Fitting in Japan and Italy

Abstract

COVID-19, which broke out globally in 2019, is an infectious disease caused by a novel strain of coronavirus, and its spread is highly contagious and concealed. Environmental vectors play an important role in viral infection and transmission, which brings new difficulties and challenges to disease prevention and control. In this paper, a type of differential equation model is constructed according to the spreading functions and characteristics of exposed individuals and environmental vectors during the virus infection process. In the proposed model, five compartments were considered, namely, susceptible individuals, exposed individuals, infected individuals, recovered individuals, and environmental vectors (contaminated with free virus particles). In particular, the re-positive factor was taken into account (i.e., recovered individuals who have lost sufficient immune protection may still return to the exposed class). With the basic reproduction number ${\mathcal{R}}_0$ of the model, the global stability of the disease-free equilibrium and uniform persistence of the model were completely analyzed. Furthermore, sufficient conditions for the global stability of the endemic equilibrium of the model were also given. Finally, the effective predictability of the model was tested by fitting COVID-19 data from Japan and Italy.

1. Introduction

In 2019, an unexplained novel coronavirus disease (COVID-19) suddenly broke out, which triggered an unprecedented public health crisis in the world. On 11 March 2020, the World Health Organization (WHO) declared COVID-19 as a global pandemic. COVID-19 is caused by a novel coronavirus named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1]. SARS-CoV-2 is regarded as the third zoonotic coronavirus emerging in the current century, after SARS-CoV in 2002 and the Middle East respiratory syndrome coronavirus (MERS-CoV) in 2012 [2,3]. Common signs of COVID-19 infection include respiratory symptoms, fever, cough, shortness of breath, and difficulty of breathing, and in more serious cases, it can lead to pneumonia, severe acute respiratory syndrome, renal failure, and even death due to alveolar damage, posing a serious threat to human life [4,5]. Among the low-risk people who died of COVID-19, 70% of them had one or more organ dysfunctions within 4 months of the initial symptoms of COVID-19, and experienced organ damage, fatigue, muscle weakness, difficulty of sleeping, anxiety, or depression after the acute phase [6,7]. On a global scale, as of 24 May 2022, countries around the world have reported more than 500 million confirmed cases of COVID-19 to the WHO, including 6.2 million deaths cases [8].

It is encouraging that in the past two years, with the joint efforts of scientists from many countries, a series of important research achievements have been made in the field of drugs for COVID-19 treatment, improving the effectiveness of vaccines, which have greatly contributed to protecting people’s health and lives [9,10].

In particular, differential equations and statistics have been widely used to construct various types of compartmental models and make rapid and effective predictions about infection trends, epidemic patterns, and key factors in the spread of infectious diseases. In addition, mathematical models based on differential equations provide a theoretical basis and optimal strategies for the prevention and treatment of infectious diseases [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27].

It is known that environmental vectors play an important role in virus infection and transmission, and bring new difficulties and challenges to the prevention and control of the disease. For example, a study reviewing 22 types of coronaviruses revealed that viruses such as SARS-CoV, SARS-CoV, MERS-CoV, and endemic human coronaviruses can persist for up to 9 days on inanimate surfaces such as metal, glass, or plastic [28]. Another experimental study found that SARS-CoV-2 can be detected on aerosols for up to 3 h, on copper for up to 4 h, on cardboard for up to 24 h, and on plastic and stainless steel for up to 3 days [29]. Moreover, SARS-CoV-2 can also be detected in domestic wastewater [30,31]. In recent years, many scholars have conducted extensive research on infectious disease models with media effects [32,33,34,35,36,37,38,39,40]. In addition, many scholars have considered the possible factor of re-positives and designed a new model based on the traditional SEIR-type to analyze the spread of the epidemic [41,42], which can also be seen in previous studies [43,44,45,46,47,48,49,50]. According to previous studies, re-positives and environmental vectors significantly impact COVID-19 dynamics and evolution, and need to be investigated widely in order to uncover the strength of disease severity and infectiousness. It has not yet been extensively explored using mathematical modeling studies. Therefore, it is necessary to use differential equation modeling methods to help understand the mechanism of environmental vectors and re-positive factors in virus transmission.

The authors of [41] consider the situation of re-infection, where recovered people may return to being considered infected people (including symptomatic and asymptomatic infected people), and these infected people can release virus particles into the environment. In order to simplify the model, we let exposed people to include lurkers and asymptomatic infected people. We assume that asymptomatic infected people are infectious, and that all infected people must undergo the incubation period.

In addition to the factors mentioned in [41], the following three types of factors are also considered: the first is that exposed individuals are still contagious, and exposed individuals can transmit viruses deposited on the surface of materials to healthy individuals through inanimate substances (public lift buttons, mail packages, etc.); the second is that the birth population is not zero, the mortality rate of different populations is also different, and recovered individuals may return to becoming susceptible individuals due to immune loss; the third is that we consider the recovered of re-positives to return to becoming exposed people.

Let the variables S(t), E(t), I(t), and R(t) represent the numbers of susceptible individuals, exposed individuals, infected individuals, and recovered individuals at time t, respectively, and M(t) represent the number of free virus particles accumulated in environmental vectors. The block diagram of the interactions among S, E, I, R, and M is shown below (Figure 1).

According to the above block diagram of the interaction between these state variables, we have the following five-dimensional ordinary differential equation model:

$$\left\{\begin{array}{cc}\hfill \dot{S}\left(t\right)=& \Lambda -{\beta }_{1}S\left(t\right)I\left(t\right)-{\beta }_{2}S\left(t\right)E\left(t\right)-{\gamma }_{3}S\left(t\right)M\left(t\right)+{\mu }_{4}R\left(t\right)-dS\left(t\right),\hfill \\ \hfill \dot{E}\left(t\right)=& {\beta }_{1}S\left(t\right)I\left(t\right)+{\beta }_{2}S\left(t\right)E\left(t\right)+{\gamma }_{3}S\left(t\right)M\left(t\right)+{\mu }_{3}R\left(t\right)-\left(d+d_1+{\mu }_1\right)E\left(t\right),\hfill \\ \hfill \dot{I}\left(t\right)=& {\mu }_{1}E\left(t\right)-\left(d+d_2+{\mu }_2\right)I\left(t\right),\hfill \\ \hfill \dot{R}\left(t\right)=& {\mu }_{2}I\left(t\right)-\left(d+{\mu }_3+{\mu }_4\right)R\left(t\right),\hfill \\ \hfill \dot{M}\left(t\right)=& {\gamma }_{1}E\left(t\right)+{\gamma }_{2}I\left(t\right)-\rho M\left(t\right).\hfill \end{array}\right.$$

(1)

The biological meanings of each parameter in model (1) are shown in

Table 1. Biological meanings of the parameters in model (1.1).

Parameter	Description	Unit
${\beta }_1$	Infection rate of infected individuals	/10,000 persons × day
${\beta }_2$	Infection rate of exposed individuals	/10,000 persons × day
${\gamma }_1$	Rate of virus release into the environment by exposed individuals	10,000 viruses/10,000 persons × day
${\gamma }_2$	Rate of virus release into the environment by infected individuals	10,000 viruses/10,000 persons × day
${\gamma }_3$	Infection rate due to environmental vectors	/10,000 viruses × day
$\rho$	Rate of decay or clearance of free virus particles in environmental vectors	/day
$\Lambda$	Population input rate	10,000 persons/day
${\mu }_1$	Rate at which exposed individuals become infected	/day
${\mu }_2$	Recovery rate of infected people	/day
${\mu }_3$	Rate of re-positives	/day
${\mu }_4$	Rate at which recovered individuals without immune protection return to the susceptible class	/day
d	Natural death rate	/day
$d_1$	Death rate of exposed individuals caused by virus	/day
$d_2$	Death rate of infected individuals caused by virus	/day

below:

Considering the biological meaning, all the parameters in Table 1 are non-negative constants, $\Lambda >0$, $d>0$, and $\rho >0$.

The main purpose of the paper can be divided into two parts. Firstly, the local and global stability of the equilibria and uniform persistence of model (1) are analyzed in detail by using the basic reproduction number. Secondly, model (1) is applied to data fitting of COVID-19 from Japan and Italy, which shows the effectiveness and predictability of the model. Furthermore, the paper conducts a sensitivity analysis of the parameters in model (1) using data from Italy, which allows us to identify the main factors related to disease transmission (for example, the infection rate of infected individuals (${\beta }_1$), infection rate of exposed individuals (${\beta }_2$), rate of virus release into the environment by infected individuals (${\gamma }_2$), infection rate due to environmental vectors (${\gamma }_3$), rate of decay or clearance of free virus particles in environmental vectors ($\rho$), rate at which exposed individuals become infected (${\mu }_1$), recovery rate of infected people (${\mu }_2$), rate of re-positives (${\mu }_3$), and the rate at which recovered individuals without immune protection return to the susceptible class (${\mu }_4$)).

2. Stability and Uniform Persistence

The initial condition of model (1) is

$$S\left(0\right)=S^0,E\left(0\right)=E^0,I\left(0\right)=I^0,R\left(0\right)=R^0,M\left(0\right)=M^0,$$

(2)

where $S^0,$$E^0,$$I^0,$$R^0$, and $M^0$ are non-negative constants.

It is easy to prove that the solution $\left(S\right(t),E(t),I(t),R(t),M(t\left)\right)$ of model (1) with the initial condition (2) is existent, unique, non-negative, and ultimately bounded in $[0,+\infty )$, and satisfies

$$\underset{t\to \infty }{lim\; sup}(S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right))\le \Lambda d^{-1},\underset{t\to \infty }{lim\; sup}M\left(t\right)\le \left({\gamma }_1+{\gamma }_2\right)\Lambda {\left(d\rho \right)}^{-1}.$$

2.1. Basic Reproduction Number and Classification of Equilibria

Model (1) always has the disease-free equilibrium $Q_0=\left(S_0,0,0,0,0\right)$, where $S_0=\Lambda d^{-1}$. By the method of the next generation matrix [51,52], we obtain the basic reproduction number of model (1), as follows:

$${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\Lambda \left(d+{\mu }_3+{\mu }_4\right)\left[\left(\rho {\beta }_2+{\gamma }_1{\gamma }_3\right)\left(d+d_2+{\mu }_2\right)+{\mu }_1\left(\rho {\beta }_1+{\gamma }_2{\gamma }_3\right)\right]}{d\rho \left[\left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-{\mu }_1{\mu }_2{\mu }_3\right]}.$$

where

$$\begin{array}{cc}& F=\left[\begin{array}{cccc}\frac{{\beta }_2\Lambda }{d}& \frac{{\beta }_1\Lambda }{d}& 0& \frac{{\gamma }_3\Lambda }{d}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{cccc}d+d_1+{\mu }_1& 0& -{\mu }_3& 0\\ -{\mu }_1& d+d_2+{\mu }_2& 0& 0\\ 0& -{\mu }_2& d+{\mu }_3+{\mu }_4& 0\\ -{\gamma }_1& -{\gamma }_2& 0& \rho \end{array}\right].\hfill \end{array}$$

It is not difficult to see that ${\mathcal{R}}_0$ can be written as ${\mathcal{R}}_0={\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3$, where

$$\begin{array}{cc}& {\mathcal{R}}_1=\frac{{\mu }_1\left(d+{\mu }_3+{\mu }_4\right){\beta }_1{S}_0}{\left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-{\mu }_1{\mu }_2{\mu }_3},\hfill \\ & {\mathcal{R}}_2=\frac{\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right){\beta }_2{S}_0}{\left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-{\mu }_1{\mu }_2{\mu }_3},\hfill \\ & {\mathcal{R}}_3=\frac{\left[{\mu }_1{\gamma }_2+{\gamma }_1\left(d+d_2+{\mu }_2\right)\right]\left(d+{\mu }_3+{\mu }_4\right){\gamma }_3{S}_0}{\rho \left[\left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-{\mu }_1{\mu }_2{\mu }_3\right]}.\hfill \end{array}$$

According to the expressions of ${\mathcal{R}}_1$, ${\mathcal{R}}_2$, and ${\mathcal{R}}_3$, it can be clearly seen that the value of the basic reproduction number ${\mathcal{R}}_0$ directly depends on the infection rate of the infected individuals (${\beta }_1$), the infection rate of exposed individuals (${\beta }_2$), and the infection rate due to environmental vectors (${\gamma }_3$).

Assuming that $(S,E,I,R,M)$ can be any equilibrium of model (1), it has the following equations:

$$\left\{\begin{array}{c}\Lambda -{\beta }_{1}SI-{\beta }_{2}SE-{\gamma }_{3}SM+{\mu }_{4}R-dS=0,\hfill \\ {\beta }_{1}SI+{\beta }_{2}SE+{\gamma }_{3}SM+{\mu }_{3}R-\left(d+d_1+{\mu }_1\right)E=0,\hfill \\ {\mu }_{1}E-\left(d+d_2+{\mu }_2\right)I=0,\hfill \\ {\mu }_{2}I-\left(d+{\mu }_3+{\mu }_4\right)R=0,\hfill \\ {\gamma }_{1}E+{\gamma }_{2}I-\rho M=0.\hfill \end{array}\right.$$

(3)

By solving the equations (3), it is easy to obtain that when ${\mathcal{R}}_0>1$, model (1) has a unique endemic equilibrium $Q^{\ast }=(S^{\ast },E^{\ast },I^{\ast },R^{\ast },M^{\ast })$, where

$$\begin{array}{cc}\hfill S^{\ast }& =\frac{\rho \left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-\rho {\mu }_1{\mu }_2{\mu }_3}{{\mu }_1\left(\rho {\beta }_1+{\gamma }_2{\gamma }_3\right)\left(d+{\mu }_3+{\mu }_4\right)+\left(\rho {\beta }_2+{\gamma }_1{\gamma }_3\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)}=\frac{\Lambda }{{\mathcal{R}}_{0}d},\hfill \\ \hfill E^{\ast }& =\frac{\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)}{{\mu }_1{\mu }_2}{R}^{\ast },I^{\ast }=\frac{\left(d+{\mu }_3+{\mu }_4\right)}{{\mu }_2}{R}^{\ast },\hfill \\ \hfill R^{\ast }& =\left(1-\frac{1}{{\mathcal{R}}_0}\right)\frac{{\mu }_1{\mu }_2\Lambda }{\left(d+d_1+{\mu }_1\right)\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)-\left({\mu }_1{\mu }_2{\mu }_3+{\mu }_1{\mu }_2{\mu }_4\right)},\hfill \\ \hfill M^{\ast }& =\frac{{\gamma }_1\left(d+d_2+{\mu }_2\right)\left(d+{\mu }_3+{\mu }_4\right)+{\mu }_1{\gamma }_2\left(d+{\mu }_3+{\mu }_4\right)}{{\mu }_1{\mu }_2\rho }{R}^{\ast }.\hfill \end{array}$$

2.2. Global Stability of the Disease-Free Equilibrium

It is clear that the set

$$\Omega =\left\{(E,I,R,M)|E,I,R,M\ge 0,S+E+I+R\le S_0,M\le \left({\gamma }_1+{\gamma }_2\right)\Lambda {\left(d\rho \right)}^{-1}\right\}$$

is attractive and positively invariant with respect to model (1). Thus, this leads to the following theorem:

Theorem 1. 

If ${\mathcal{R}}_0<1$, the disease-free equilibrium $Q_0$ of model (1) is globally asymptotically stable with respect to Ω.

Proof. 

Firstly, let us show that the disease-free equilibrium $Q_0$ is locally asymptotically stable. For convenience, let $d+d_1+{\mu }_1=a_1,d+d_2+{\mu }_2=a_2$, and $d+{\mu }_3+{\mu }_4=a_3$. Hence, $a_1{a}_2{a}_3-{\mu }_1{\mu }_2{\mu }_3>0$.

Using simple calculations, we can determine that the characteristic equation of model (1) at the disease-free equilibrium $Q_0$ is $(\lambda +d)g\left(\lambda \right)=0$, where

$$\begin{array}{cc}\hfill g\left(\lambda \right)=& -\frac{{\gamma }_3\Lambda }{d}(\lambda +a_3)[{\mu }_1{\gamma }_2+{\gamma }_1(\lambda +a_2)]+(\lambda +\rho )(\lambda -\frac{{\beta }_2\Lambda }{d}+a_1)(\lambda +a_2)(\lambda +a_3)\hfill \\ \hfill & -(\lambda +\rho )[\frac{{\mu }_1{\beta }_1\Lambda }{d}(\lambda +a_3)+{\mu }_1{\mu }_2{\mu }_3].\hfill \end{array}$$

Clearly, ${\lambda }_1=-d$ is a negative root. Let us further show that when ${\mathcal{R}}_0<1$, all the roots of the equation $g\left(\lambda \right)=0$ have negative real parts. Suppose that $\lambda =x+iy$$(x\ge 0)$ is an eigenvalue of the equation $g\left(\lambda \right)=0$.

Then, the equation $g\left(\lambda \right)=0$ can be rewritten in the following equivalent form:

$$\begin{array}{ccc}\hfill 1=& \frac{{\mu }_1{\gamma }_2{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)(\lambda +a_2)}+\frac{{\gamma }_1{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)}+\frac{{\beta }_2\Lambda }{d(\lambda +a_1)}\hfill & \hfill +\frac{{\mu }_1{\beta }_1\Lambda }{d(\lambda +a_1)(\lambda +a_2)}\\ \hfill & +\frac{{\mu }_1{\mu }_2{\mu }_3}{(\lambda +a_1)(\lambda +a_2)(\lambda +a_3)}.\hfill \end{array}$$

Taking the modulus on both sides of the above equality, it becomes

$$\begin{array}{c}1=|\frac{{\mu }_1{\gamma }_2{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)(\lambda +a_2)}+\frac{{\gamma }_1{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)}+\frac{{\beta }_2\Lambda }{d(\lambda +a_1)}+\frac{{\mu }_1{\beta }_1\Lambda }{d(\lambda +a_1)(\lambda +a_2)}\hfill \\ +\frac{{\mu }_1{\mu }_2{\mu }_3}{(\lambda +a_1)(\lambda +a_2)(\lambda +a_3)}|\hfill \\ \le |\frac{{\mu }_1{\gamma }_2{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)(\lambda +a_2)}|+|\frac{{\gamma }_1{\gamma }_3\Lambda }{d(\lambda +\rho )(\lambda +a_1)}|+|\frac{{\beta }_2\Lambda }{d(\lambda +a_1)}|+|\frac{{\mu }_1{\beta }_1\Lambda }{d(\lambda +a_1)(\lambda +a_2)}|\hfill \\ +|\frac{{\mu }_1{\mu }_2{\mu }_3}{(\lambda +a_1)(\lambda +a_2)(\lambda +a_3)}|\hfill \\ \le \frac{{\mu }_1{\gamma }_2{\gamma }_3\Lambda }{d\rho a_1{a}_2}+\frac{{\gamma }_1{\gamma }_3{a}_2\Lambda }{d\rho a_1{a}_2}+\frac{{\beta }_2\rho a_2\Lambda }{d\rho a_1{a}_2}+\frac{{\mu }_1{\beta }_1\rho \Lambda }{d\rho a_1{a}_2}+\frac{{\mu }_1{\mu }_2{\mu }_3}{a_1{a}_2{a}_3}\hfill \\ =\frac{\Lambda ({\mu }_1{\gamma }_2{\gamma }_3+{\mu }_1{\beta }_1\rho +{\gamma }_1{\gamma }_3{a}_2+{\beta }_2\rho a_2)}{d\rho a_1{a}_2}+\frac{{\mu }_1{\mu }_2{\mu }_3}{a_1{a}_2{a}_3}\hfill \\ =\frac{(a_1{a}_2{a}_3-{\mu }_1{\mu }_2{\mu }_3)}{a_1{a}_2{a}_3}{\mathcal{R}}_0+\frac{{\mu }_1{\mu }_2{\mu }_3}{a_1{a}_2{a}_3}\hfill \\ ={\mathcal{R}}_0+(1-{\mathcal{R}}_0)\frac{{\mu }_1{\mu }_2{\mu }_3}{a_1{a}_2{a}_3}.\hfill \end{array}$$

Note that $a_1{a}_2{a}_3-{\mu }_1{\mu }_2{\mu }_3>0$ and ${\mathcal{R}}_0<1$, meaning that the above inequality is not valid. This proves that when ${\mathcal{R}}_0<1$, all the roots of the equation $g\left(\lambda \right)=0$ have negative real parts. Hence, the disease-free equilibrium $Q_0$ is locally asymptotically stable.

Next, we show that the disease-free equilibrium $Q_0$ is globally attractable. In fact, for $t\ge 0$, model (1) implies

$$\left\{\begin{array}{c}\dot{E}\left(t\right)\le {\beta }_1{S}_{0}I\left(t\right)+{\beta }_2{S}_{0}E\left(t\right)+{\gamma }_3{S}_{0}M\left(t\right)+{\mu }_{3}R\left(t\right)-\left(d+d_1+{\mu }_1\right)E\left(t\right),\hfill \\ \dot{I}\left(t\right)\le {\mu }_{1}E\left(t\right)-\left(d+d_2+{\mu }_2\right)I\left(t\right),\hfill \\ \dot{R}\left(t\right)\le {\mu }_{2}I\left(t\right)-\left(d+{\mu }_3+{\mu }_4\right)R\left(t\right),\hfill \\ \dot{M}\left(t\right)\le {\gamma }_{1}E\left(t\right)+{\gamma }_{2}I\left(t\right)-\rho M\left(t\right).\hfill \end{array}\right.$$

Let $Y={(y_1,y_2,y_3,y_4)}^T$, and consider the comparison system $\frac{dY}{dt}=(F-V)Y$. As the condition ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)<1$ means that all eigenvalues of matrix $F-V$ have negative real parts, the trivial solution of the comparison system is asymptotically stable. Therefore, it follows that $E\left(t\right)\le y_1\left(t\right)\to 0$, $I\left(t\right)\le y_2\left(t\right)\to 0$, $R\left(t\right)\le y_3\left(t\right)\to 0$ and $M\left(t\right)\le y_4\left(t\right)\to 0$ for $t\to +\infty$. Furthermore, from the first equation of model (1), it is easy to obtain $S\left(t\right)\to S_0$ for $t\to +\infty$. This proves that the disease-free equilibrium $Q_0$ is globally attractive with respect to the set $\Omega$. □

2.3. Uniform Persistence and Global Stability of the Endemic Equilibrium

In the subsection, there is the following Theorem 2 for the uniform persistence of model (1).

Theorem 2. 

If ${\mathcal{R}}_0>1$, then model (1) is uniformly persistent, and each positive solution ${(S\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right),M\left(t\right))}^T$ of model (1) with the initial condition (2) satisfies

$$\begin{array}{cc}& \underset{t\to \infty }{lim\; inf}S\left(t\right)\ge \frac{\Lambda d}{\Lambda \left[{\beta }_1+{\beta }_2+{\gamma }_3\left({\gamma }_1+{\gamma }_2\right){\rho }^{-1}\right]+d^2}\equiv v_1,\underset{t\to \infty }{lim\; inf}E\left(t\right)\ge \delta E^{\ast }{e}^{-T{a}_1}\equiv v_2,\hfill \\ & \underset{t\to \infty }{lim\; inf}I\left(t\right)\ge \frac{{\mu }_1}{a_2}{v}_2,\underset{t\to \infty }{lim\; inf}R\left(t\right)\ge \frac{{\mu }_1{\mu }_2}{a_2{a}_3}{v}_2,\underset{t\to \infty }{lim\; inf}M\left(t\right)\ge \left(\frac{{\gamma }_1}{\rho }+\frac{{\gamma }_2{\mu }_1}{a_2\rho }\right)v_2,\hfill \end{array}$$

where $\delta >0$ and $T>0$, and satisfies

$$q\equiv \left(\frac{{\beta }_1{v}_2}{a_2}+{\beta }_2+\frac{{\gamma }_1{\gamma }_3}{\rho }+\frac{{\gamma }_2{\gamma }_3{\mu }_1}{a_2\rho }\right)\delta E^{\ast }+d,\frac{\Lambda }{q}>S^{\ast },S^{\Delta }\equiv \frac{\Lambda }{q}\left(1-{\mathrm{e}}^{-Tq}\right)>S^{\ast }.$$

To complete the proof of Theorem 2, the key point is to show that the estimation $\underset{t\to +\infty }{lim\; inf}E\left(t\right)\ge v_2$ holds. The detailed proof is similar to [53,54,55,56] and has been omitted here.

From a biological point of view, Theorem 2 indicates that as long as the basic reproduction number ${\mathcal{R}}_0>1$, the disease infection cannot be eliminated and will permanently exist.

Next, let us consider the global asymptotic stability of the endemic equilibrium $Q^{\ast }$ of model (1). To simplify model (1), it is assumed that the death rate of the exposed and infected individuals caused by the virus is zero, i.e., $\left(H1\right)$$d_1=d_2=0$.

Let $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$, then for $t\ge 0$, we have $\dot{N}\left(t\right)=\Lambda -dN\left(t\right)$, which implies $\underset{t\to +\infty }{lim}N\left(t\right)=S_0$. On the hyperplane $N=S_0$, model (1) can be transformed into the following equivalent four-dimensional system:

$$\left\{\begin{array}{cc}\hfill \dot{S}\left(t\right)=& \Lambda -{\beta }_{1}S\left(t\right)I\left(t\right)-{\beta }_{2}S\left(t\right)E\left(t\right)-{\gamma }_{3}S\left(t\right)M\left(t\right)+{\mu }_4(S_0-S\left(t\right)-E\left(t\right)-I\left(t\right))-dS\left(t\right),\hfill \\ \hfill \dot{E}\left(t\right)=& {\beta }_{1}S\left(t\right)I\left(t\right)+{\beta }_{2}S\left(t\right)E\left(t\right)+{\gamma }_{3}S\left(t\right)M\left(t\right)+{\mu }_3(S_0-S\left(t\right)-E\left(t\right)-I\left(t\right))-\left(d+{\mu }_1\right)E\left(t\right),\hfill \\ \hfill \dot{I}\left(t\right)=& {\mu }_{1}E\left(t\right)-\left(d+{\mu }_2\right)I\left(t\right),\hfill \\ \hfill \dot{M}\left(t\right)=& {\gamma }_{1}E\left(t\right)+{\gamma }_{2}I\left(t\right)-\rho M\left(t\right).\hfill \end{array}\right.$$

(4)

In addition, the following conditions are also used:

$\left(H2\right)$${\mu }_2+{\mu }_3-2{\mu }_1-{\mu }_4>0,d+{\mu }_1+{\mu }_3-\left(2{\beta }_1+2{\beta }_2\right)S_0-\rho -2{\mu }_4>0,d+{\mu }_2-{\mu }_1>0$, ${\gamma }_1\ge {\gamma }_2$, $2{\mu }_4\ge {\mu }_1+{\mu }_3.$

By using a similar method to what was used in [57,58,59], we have Theorem 3.

Theorem 3. 

If the conditions $\left(H1\right)$–$\left(H2\right)$ hold, then, when ${\mathcal{R}}_0>1$, the endemic equilibrium $Q^{\ast }$ of system (4) is globally asymptotically stable with respect to Ω.

Proof. 

For system (4), the Jacobian matrix J and second additive compound matrix $J^{\left[2\right]}$ are given by

$$J=\left(\begin{array}{cccc}-{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-d-{\mu }_4& -{\beta }_{2}S-{\mu }_4& -{\beta }_{1}S-{\mu }_4& -{\gamma }_{3}S\\ {\beta }_{1}I+{\beta }_{2}E+{\gamma }_{3}M-{\mu }_3& {\beta }_{2}S-d-{\mu }_1-{\mu }_3& {\beta }_{1}S-{\mu }_3& {\gamma }_{3}S\\ 0& {\mu }_1& -d-{\mu }_2& 0\\ 0& {\gamma }_1& {\gamma }_2& -\rho \end{array}\right)$$

and

$$J^{\left[2\right]}=\left(\begin{array}{cccccc}{M}_{11}& {\beta }_{1}S-{\mu }_3& {\gamma }_{3}S& {\beta }_{1}S+{\mu }_4& {\gamma }_{3}S& 0\\ {\mu }_1& M_{22}& 0& -{\beta }_{2}S-{\mu }_4& 0& {\gamma }_{3}S\\ {\gamma }_1& {\gamma }_2& M_{33}& 0& -{\beta }_{2}S-{\mu }_4& -{\beta }_{1}S-{\mu }_4\\ 0& q& 0& M_{44}& 0& -{\gamma }_{3}S\\ 0& 0& q& {\gamma }_2& M_{55}& {\beta }_{1}S-{\mu }_3\\ 0& 0& 0& -{\gamma }_1& {\mu }_1& M_{66}\end{array}\right),$$

where

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{11}=-{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-d-{\mu }_4+{\beta }_{2}S-d-{\mu }_1-{\mu }_3$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{22}=-{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-d-{\mu }_4-d-{\mu }_2$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{33}=-{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-d-{\mu }_4-\rho$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{44}={\beta }_{2}S-2d-{\mu }_1-{\mu }_3-{\mu }_2$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{55}={\beta }_{2}S-d-{\mu }_1-{\mu }_3-\rho$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{M}_{66}=-d-{\mu }_2-\rho$,

$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}q={\beta }_{1}I+{\beta }_{2}E+{\gamma }_{3}M-{\mu }_3$.

Let

$$P=P(S,E,I,M)=\left(\begin{array}{cccccc}\frac{1}{E}& 0& 0& 0& 0& 0\\ 0& \frac{1}{E}& 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{E}& 0& 0\\ 0& 0& \frac{1}{M}& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{M}& 0\\ 0& 0& 0& 0& 0& \frac{1}{M}\end{array}\right),$$

using f to represent the vector form of system (4). Then, $P_f$ is the directional derivative of $P\left(x\right)$ along the direction of f. It follows that $P_f{P}^{-1}=diag\left\{-\frac{\dot{E}}{E},-\frac{\dot{E}}{E},-\frac{\dot{E}}{E},-\frac{\dot{M}}{M},-\frac{\dot{M}}{M},-\frac{\dot{M}}{M}\right\}$, and

$$\begin{array}{cc}\hfill Q(S,E,I,M)=& P_f{P}^{-1}+P{J}^{\left[2\right]}{P}^{-1}\hfill \\ \hfill =& \left(\begin{array}{cccccc}{M}_{11}-\frac{\dot{E}}{E}& {\beta }_{1}S-{\mu }_3& {\beta }_{1}S+{\mu }_4& \frac{{\gamma }_{3}SM}{E}& \frac{{\gamma }_{3}SM}{E}& 0\\ {\mu }_1& M_{22}-\frac{\dot{E}}{E}& -{\beta }_{2}S-{\mu }_4& 0& 0& \frac{{\gamma }_{3}SM}{E}\\ 0& q& M_{44}-\frac{\dot{E}}{E}& 0& 0& -\frac{{\gamma }_{3}SM}{E}\\ \frac{{\gamma }_{2}E}{M}& \frac{{\gamma }_{2}E}{M}& 0& M_{33}-\frac{\dot{M}}{M}& -{\beta }_{2}S-{\mu }_4& -{\beta }_{1}S-{\mu }_4\\ 0& 0& \frac{{\gamma }_{2}E}{M}& q& M_{55}-\frac{\dot{M}}{M}& {\beta }_{1}S-{\mu }_3\\ 0& 0& -\frac{{\gamma }_{1}E}{M}& 0& {\mu }_1& M_{55}-\frac{\dot{M}}{M}\end{array}\right).\hfill \end{array}$$

The matrix $Q\left(S,E,I,M\right)$ can be written in block form: $Q\left(S,E,I,M\right)={\left(Q_{ij}\right)}_{4\times 4}$, where

$$\begin{array}{cc}\hfill & Q_{11}=M_{11}-\frac{\dot{E}}{E},Q_{12}=\left({\beta }_{1}S-{\mu }_3,{\beta }_{1}S+{\mu }_4\right),Q_{13}=\left(\frac{{\gamma }_{3}SM}{E},\frac{{\gamma }_{3}SM}{E}\right),Q_{14}=0,\hfill \\ \hfill & Q_{21}={\left({\mu }_1,0\right)}^T,Q_{22}=\left(\begin{array}{cc}{M}_{22}-\frac{\dot{E}}{E}& -{\beta }_{2}S-{\mu }_4\\ q& M_{44}-\frac{\dot{E}}{E}\end{array}\right),Q_{23}=\left(\begin{array}{cc}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \end{array}\right),\hfill \\ \hfill & Q_{24}={\left(\frac{{\gamma }_{3}SM}{E},-\frac{{\gamma }_{3}SM}{E}\right)}^T,Q_{31}={\left(\frac{{\gamma }_{2}E}{M},0\right)}^T,Q_{32}=\left(\begin{array}{cc}\frac{{\gamma }_{2}E}{M}& 0\\ 0& \frac{{\gamma }_{2}E}{M}\end{array}\right),\hfill \\ \hfill & Q_{33}=\left(\begin{array}{cc}{M}_{33}-\frac{\dot{M}}{M}& -{\beta }_{2}S-{\mu }_4\\ q& M_{55}-\frac{\dot{M}}{M}\end{array}\right),Q_{34}={\left(-{\beta }_{1}S-{\mu }_4,{\beta }_{1}S-{\mu }_3\right)}^T,Q_{41}=0,\hfill \\ \hfill & Q_{42}=(0,-\frac{{\gamma }_{1}E}{M}),Q_{43}=(0,{\mu }_1),Q_{44}=M_{66}-\frac{\dot{M}}{M}.\hfill \end{array}$$

Hence, we have the following estimation (for an example, see [57,58,59]),

$$\sigma \left(Q(S,E,I,M)\right)\le sup\left\{g_1,g_2,g_3,g_4\right\},$$

where $g_1={\sigma }_1\left(Q_{11}\right)+\left|Q_{12}\right|+\left|Q_{13}\right|+\left|Q_{14}\right|,$ $g_2={\sigma }_1\left(Q_{22}\right)+\left|Q_{21}\right|+\left|Q_{23}\right|+\left|Q_{24}\right|,$ $g_3={\sigma }_1\left(Q_{33}\right)+\left|Q_{31}\right|+\left|Q_{32}\right|+\left|Q_{34}\right|,$$g_4={\sigma }_1\left(Q_{44}\right)+\left|Q_{41}\right|+\left|Q_{42}\right|+\left|Q_{43}\right|.$ $\left|Q_{ij}\right|(i\ne j,i,j=1,2,3,4)$ are matrix norms with respect to the $l_1$ vector norm, and ${\sigma }_1$ denotes the Lozinski $\breve{ı}$ measures with respect to the $l_1$ norm. Furthermore, it is easy to obtain that

$$\begin{array}{cc}\hfill \sigma \left(Q_{11}\right)=& -{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-2d-{\mu }_4+{\beta }_{2}S-{\mu }_1-{\mu }_3-\frac{\dot{E}}{E},\hfill \\ \hfill \sigma \left(Q_{22}\right)\le & -{\mu }_2-2d-\frac{\dot{E}}{E}+max\left\{-{\mu }_4,2{\beta }_{2}S+{\mu }_4-{\mu }_1-{\mu }_3\right\},\hfill \\ \hfill \sigma \left(Q_{33}\right)\le & -\rho -d-\frac{\dot{M}}{M}+max\left\{-{\mu }_4,2{\beta }_{2}S+{\mu }_4-{\mu }_1-{\mu }_3\right\},\hfill \\ \hfill \sigma \left(Q_{44}\right)=& -d-{\mu }_2-\rho -\frac{\dot{M}}{M},\left|Q_{12}\right|={\beta }_{1}S+{\mu }_4,\left|Q_{13}\right|=\frac{{\gamma }_{3}SM}{E},\hfill \end{array}$$

$$\left|Q_{14}\right|=0,\left|Q_{21}\right|={\mu }_1,\left|Q_{23}\right|=0,\left|Q_{24}\right|=\frac{2{\gamma }_{3}SM}{E},\left|Q_{31}\right|=\frac{{\gamma }_{2}E}{M},$$

$$\left|Q_{32}\right|=\frac{{\gamma }_{2}E}{M},\left|Q_{34}\right|<2{\beta }_{1}S+{\mu }_4,\left|Q_{41}\right|=0,\left|Q_{42}\right|=\frac{{\gamma }_{1}E}{M},\left|Q_{43}\right|={\mu }_1.$$

From condition $\left(H2\right)$ and system (4), we have

$$max\left\{-{\mu }_4,2{\beta }_{2}S+{\mu }_4-{\mu }_1-{\mu }_3\right\}=2{\beta }_{2}S+{\mu }_4-{\mu }_1-{\mu }_3,$$

and

$$\begin{array}{cc}\hfill \frac{\dot{E}}{E}=& \frac{{\beta }_{1}SI}{E}+{\beta }_{2}S+\frac{{\gamma }_{3}SM}{E}+\frac{{\mu }_{3}R}{E}-(d+{\mu }_1)>\frac{{\beta }_{1}SI}{E}+{\beta }_{2}S+\frac{{\gamma }_{3}SM}{E}-(d+{\mu }_1),\hfill \\ \hfill \frac{\dot{M}}{M}=& \frac{{\gamma }_{1}E}{M}+\frac{{\gamma }_{2}I}{M}-\rho >\frac{{\gamma }_{1}E}{M}-\rho \ge \frac{{\gamma }_{2}E}{M}-\rho .\hfill \end{array}$$

Hence, we have

$$\begin{array}{cc}\hfill g_1=& -{\beta }_{1}I-{\beta }_{2}E-{\gamma }_{3}M-2d-{\mu }_4+{\beta }_{2}S-{\mu }_1-{\mu }_3-\frac{\dot{E}}{E}+{\beta }_{1}S+{\mu }_4+\frac{{\gamma }_{3}SM}{E}\hfill \\ \hfill <& -{\beta }_{1}I-{\beta }_{2}S-{\gamma }_{3}M-\frac{{\beta }_{1}SI}{E}-d-{\mu }_3+{\beta }_{1}S\hfill \\ \hfill <& \frac{{\beta }_1\Lambda }{d}-d-{\mu }_3\equiv -{\theta }_1,\hfill \\ \hfill g_2\le & -{\mu }_2-2d-{\mu }_3-\frac{\dot{E}}{E}+2{\beta }_{2}S+{\mu }_4-{\mu }_1+{\mu }_1+\frac{2{\gamma }_{3}SM}{E}\hfill \\ \hfill <& -{\mu }_2-d-{\mu }_3+\frac{{\gamma }_{3}SM}{E}+{\beta }_{2}S-\frac{{\beta }_{1}SI}{E}+{\mu }_1+{\mu }_4\hfill \\ \hfill <& \frac{\dot{E}}{E}+2{\mu }_1+{\mu }_4-{\mu }_2-{\mu }_3\equiv \frac{\dot{E}}{E}-{\theta }_2,\hfill \\ \hfill g_3<& -\rho -{\mu }_3-d-\frac{\dot{M}}{M}+2{\beta }_{2}S+{\mu }_4-{\mu }_1+\frac{2{\gamma }_{2}E}{M}+2{\beta }_{1}S+{\mu }_4\hfill \\ \hfill <& \frac{{\gamma }_{2}E}{M}-\rho +\rho -{\mu }_3+\left(2{\beta }_1+2{\beta }_2\right)S+2{\mu }_4-d-{\mu }_1\hfill \\ \hfill <& \frac{\dot{M}}{M}+\left(2{\beta }_1+2{\beta }_2\right)\frac{\Lambda }{d}+\rho +2{\mu }_4-d-{\mu }_1-{\mu }_3\equiv \frac{\dot{M}}{M}-{\theta }_3,\hfill \\ \hfill g_4=& -d-{\mu }_2-\rho -\frac{\dot{M}}{M}+\frac{{\gamma }_{1}E}{M}+{\mu }_1\hfill \\ \hfill <& {\mu }_1-d-{\mu }_2\equiv -{\theta }_4.\hfill \end{array}$$

Again, from condition $\left(H2\right)$, we can obtain

$${\theta }_3<{\theta }_1,\overline{b}=min\left\{{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4\right\}=min\left\{{\theta }_2,{\theta }_3,{\theta }_4\right\}>0,$$

and

$$g_1\le -\overline{b},g_2<\frac{\dot{E}}{E}-\overline{b},g_3<\frac{\dot{M}}{M}-\overline{b},g_4\le -\overline{b}.$$

Note that when ${\mathcal{R}}_0>1$, system (4) is uniformly persistent. Using a similar argument as in [59], we find that the endemic equilibrium $Q^{\ast }$ of system (4) is globally asymptotically stable. □

Take $\Lambda =1000$, ${\beta }_1=5\times {10}^{-7}$, ${\beta }_2=5\times {10}^{-7}$, ${\gamma }_1=2\times {10}^{-3}$, ${\gamma }_2=2\times {10}^{-3}$, ${\gamma }_3=2\times {10}^{-3}$, ${\mu }_1=0.2$, ${\mu }_2=0.5$, ${\mu }_3=0.1$, ${\mu }_4=0.15$, $d=0.1$, $d_1=0.1,d_2=0.1,$ and $\rho =0.05$. Using numerical calculations, we obtain the basic reproduction number $R_0=2.625>1$, and the inequalities in $\left(H1\right)$–$\left(H2\right)$ hold.

3. Applications of the Model in Japan and Italy

In this section, we use model (1) to fit the data of confirmed cases and recovered cases of COVID-19 in Japan and Italy (data were taken from the Johns Hopkins University Center for Systems Science and Engineering, https://github.com/CSSEGISandData/COVID-19 accessed on 1 July 2022), and make short-term predictions about disease infection trends. Meanwhile, based on the basic reproduction number $R_0$ and data in Italy, we carried out a sensitivity analysis and the main factors related to disease transmission were captured.

3.1. Predicted Cases for Cumulative Confirmed and Recovered Cases Based on Data in Japan and Italy

First, from [3,60,61,62] and data published by the World Health Organization (WHO), the parameters $\frac{1}{{\mu }_1}$ (incubation period), $\frac{1}{{\mu }_2}$ (infectious period), $\frac{1}{{\mu }_4}$ (time of immune protection), $\frac{1}{d}$ (mean lifetime), and $d_2$ (death rate of the infected individuals caused by the virus) in model (1) can take the values as shown in

Table 2. The estimated values of parameters in the model.

Para-Meter	Definitions	Value (Japan)	Value (Italy)	Unit	Source
${\beta }_1$	See Table 1	$2.744\times {10}^{-7}$	$3.195\times {10}^{-6}$	/10,000 persons × day	LSM
${\beta }_2$	See Table 1	$6.766\times {10}^{-7}$	$5.123\times {10}^{-6}$	/10,000 persons × day	LSM
${\gamma }_1$	See Table 1	0.6513	0.3966	10,000 viruses/10,000 persons × day	LSM
${\gamma }_2$	See Table 1	0.2791	0.6091	10,000 viruses/10,000 persons × day	LSM
${\gamma }_3$	See Table 1	$1.559\times {10}^{-8}$	$5.235\times {10}^{-7}$	/10,000 viruses × day	LSM
$\rho$	See Table 1	0.8901	0.5957	/day	LSM
$\Lambda$	See Table 1	0.3916	0.3358	10,000 persons/day	LSM
$\frac{1}{{\mu }_1}$	Incubation period	5.2	5.2	day	[3]
$\frac{1}{{\mu }_2}$	Infectious period	20	20	day	[60]
${\mu }_3$	See Table 1	0.03051	0.09186	/day	LSM
$\frac{1}{{\mu }_4}$	Time of immune protection	180	180	day	[61]
$\frac{1}{d}$	Mean lifetime	$84\times 365$	$83\times 365$	day	WHO
$d_1$	See Table 1	$2.025\times {10}^{-4}$	$3.371\times {10}^{-4}$	/day	LSM
$d_2$	See Table 1	$1.425\times {10}^{-3}$	$2.512\times {10}^{-3}$	/day	[62]

below. Furthermore, by using data of COVID-19 in Japan (20 May–18 June 2022) and Italy (12 January–10 February 2022) and the least squares method (LSM), the remaining parameters in model (1) take values as shown in Table 2 below.

Based on the parameter values given in Table 2, we obtain Figure 2a,b, which shows that model (1) fits well with the evolution of cumulative confirmed cases and recovered cases in Japan from 20 May–18 June 2022. Similarly, we have Figure 3a,b that shows that model (1) also fits well with the evolution of the cumulative confirmed cases and recovered cases in Italy from 12 January–10 February 2022.

Next, we use model (1) and the parameter values given in Table 2 to predict the cumulative confirmed cases and recovered cases in Japan (from 19–28 June 2022) and Italy (from 11–20 February 2022). These are shown in

Table 3. Predicted values and reported data for confirmed cases in Japan.

Data	Cumulative Confirmed Cases (Japan)
	Reported	Predicted
$6.19$	$917.2745$	$912.8151$
$6.20$	$920.0541$	$913.6510$
$6.21$	$923.5918$	$914.4373$
$6.22$	$927.3191$	$915.1782$
$6.23$	$930.9861$	$915.8729$
$6.24$	$931.5674$	$916.5263$
$6.25$	$933.2259$	$917.1391$
$6.26$	$934.6495$	$917.7133$
$6.27$	$935.6060$	$918.2511$
$6.28$	$936.5433$	$918.7538$

,

Table 4. Predicted values and reported data for recovered cases in Japan.

Data	Cumulative Recovered Cases (Japan)
	Reported	Predicted
$6.19$	$897.7861$	$894.7521$
$6.20$	$900.3232$	$895.9912$
$6.21$	$903.6874$	$897.1613$
$6.22$	$905.2890$	$898.2644$
$6.23$	$908.7801$	$899.3047$
$6.24$	$910.2259$	$900.2809$
$6.25$	$914.0274$	$901.1989$
$6.26$	$916.9366$	$902.0616$
$6.27$	$917.9099$	$902.8693$
$6.28$	$918.6873$	$903.0266$

,

Table 5. Predicted values and reported data for confirmed cases in Italy.

Data	Cumulative Confirmed Cases (Italy)
	Reported	Predicted
$2.11$	$1201.3631$	$1191.4833$
$2.12$	$1212.1109$	$1198.4161$
$2.13$	$1220.3330$	$1203.0440$
$2.14$	$1228.5675$	$1209.3611$
$2.15$	$1236.4451$	$1214.3588$
$2.16$	$1242.5474$	$1218.0709$
$2.17$	$1250.5343$	$1223.4792$
$2.18$	$1255.3398$	$1226.6014$
$2.19$	$1260.7098$	$1230.4304$
$2.20$	$1266.1773$	$1233.9731$

and

Table 6. Predicted values and reported data for recovered cases in Italy.

Data	Cumulative Recovered Cases (Italy)
	Reported	Predicted
$2.11$	$908.0136$	$904.2121$
$2.12$	$918.9429$	$911.0671$
$2.13$	$927.6892$	$917.8210$
$2.14$	$934.5987$	$923.4734$
$2.15$	$940.2540$	$928.2089$
$2.16$	$946.3380$	$933.3926$
$2.17$	$953.3268$	$938.6280$
$2.18$	$960.2908$	$942.6809$
$2.19$	$966.0380$	$946.5611$
$2.20$	$970.5147$	$949.2560$

.

From Table 3, Table 4, Table 5 and Table 6, it can be seen that the predicted values for confirmed and recovered cases are within a range of 97 to $103\%$ of the reported data. With the increase in time t, the relative error also increases.

Furthermore, we use MAPE (mean absolute percentage error) and RMSPE (root mean square error) to assess the reliability of model (1) in data fitting (for examples, see [63,64]). For convenience, we introduce the definitions of MAPE and RMSPE and their evaluation standard, as follows (

Table 7. MAPE/RMSPE evaluation standard.

MAPE/RMSPE	Predictive Ability
<10%	Precision prediction
10–20%	Good prediction
20–50%	Reasonable prediction
>50%	Unreasonable prediction

):

$$\mathrm{MAPE}=\frac{{\displaystyle \sum _{i=1}^n}\left|\frac{y_i-y_i^{\prime }}{y_i}\right|}{n}\times 100\%,\mathrm{RMSPE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n}{\left[\frac{y_i-y_i^{\prime }}{y_i}\right]}^2}{n-1}}\times 100\%,$$

Here, $y_i$ and $y_i^{\prime }$ are the reported data and predicted data, respectively, and the positive integer n is the number of predicted data. Using the data in Table 3, Table 4, Table 5 and Table 6, we obtain

Table 8. MAPE/RMSPE values and their predictive ability.

Data Type	MAPE	Predictive Ability	RMSPE	Predictive Ability
Confirmed cases in Japan	$1.40\%$	Precision prediction	$1.56\%$	Precision prediction
Reovered cases in Japan	$1.08\%$	Precision prediction	$1.24\%$	Precision prediction
Confirmed cases in Italy	$1.81\%$	Precision prediction	$1.99\%$	Precision prediction
Recovered cases in Italy	$1.38\%$	Precision prediction	$1.55\%$	Precision prediction

.

3.2. Sensitivity Analysis Based on the Basic Reproduction Number and Data from Italy

From Theorem 1 and 2 we see that if the basic reproduction number ${\mathcal{R}}_0<1$ (or ${\mathcal{R}}_0>1$), the disease infection can be cleared (or will become an endemic disease). Therefore, it is very necessary to make a sensitivity analysis based on the basic reproduction number ${\mathcal{R}}_0$, and then capture the main factors related to disease transmission.

The normalized sensitivity index of ${\mathcal{R}}_0$ is defined as ${\gamma }_{ϵ}^{{\mathcal{R}}_0}\equiv \frac{\partial {\mathcal{R}}_0}{\partial ϵ}\frac{ϵ}{{\mathcal{R}}_0}.$ Here, $ϵ$ can be regarded as any parameter in ${\mathcal{R}}_0$ [65]. By using Italy’s parameter values in Table 3, we obtain

Table 9. The normalized sensitivity index based on ${\mathcal{R}}_0$ and data from Italy.

Parameter	${\beta }_1$	${\beta }_2$	${\mu }_1$	${\mu }_2$	${\mu }_3$	${\mu }_4$
Sensitivity Index	0.5827	0.2993	−0.2501	−0.1816	0.2618	−0.0833
Parameter	${\gamma }_1$	${\gamma }_2$	${\gamma }_3$	$\rho$	$d_1$	$d_2$
Sensitivity Index	0.0200	0.0979	0.1180	−0.1180	−0.0041	−0.0713

and Figure 4.

From Table 9 and Figure 4, we first observe that the parameters ${\beta }_1$ (the infection rate of infected individuals) and ${\beta }_2$ (the infection rate of exposed individuals) are the most sensitive (and positively correlated) parameters to the basic reproduction number ${\mathcal{R}}_0$. When ${\beta }_1$ increases by $10\%$ while other parameters remain unchanged, $R_0$ increases by $5.827\%$. Similarly, when ${\beta }_2$ increases by $10\%$ while other parameters remain unchanged, $R_0$ increases by $2.993\%$. Moreover, the parameters ${\mu }_1$ (the rate at which exposed individuals become infected) and t${\mu }_2$ (the recovery rate of infected people) also have high sensitivities (but are negatively correlated) to the basic reproduction number ${\mathcal{R}}_0$. When ${\mu }_1$ increases by $10\%$ while other parameters remain unchanged, $R_0$ decreases by $2.501\%$. Similarly, when ${\mu }_2$ increases by $10\%$ while other parameters remain unchanged, $R_0$ decreases by $1.816\%$.

On the other hand, we also see from Table 9 and Figure 4 that the basic reproduction number ${\mathcal{R}}_0$ has a strong dependence (positively/negatively correlated) on environmental and re-positive factors. Specifically, the parameters ${\mu }_3$ (the rate of re-positives), ${\gamma }_2$ (the rate of virus release into the environment by infected individuals), ${\gamma }_3$ (the infection rate due to environmental vectors), and $\rho$ (the rate at which the free virus particles decay or are cleared in environmental vectors) are also more sensitive to the basic reproduction number ${\mathcal{R}}_0$. When ${\gamma }_2$, ${\gamma }_3$, and ${\mu }_3$ increases by $10\%$ while other parameters remain unchanged, $R_0$ increases by 0.979, 1.180 and $2.618\%$, respectively; however, when $\rho$ increases by $10\%$, $R_0$ decreases by $1.180\%$.

This suggests that we should strengthen prevention and control to reduce the transmission rate of the disease, shorten the treatment cycle, increase antibodies by vaccination, pay attention to environmental sanitation, and reduce the risk of environmental infection.

At the end of this subsection, we present a discussion of the control strategies related to environmental and re-positive factors using the basic reproduction number ${\mathcal{R}}_0$.

Consider the basic reproduction number ${\mathcal{R}}_0$ as a function ${\mathcal{R}}_0={\mathcal{R}}_0({\mu }_3,{\gamma }_3)$ with respect to the parameters ${\mu }_3$ and ${\gamma }_3$, and let the values of all the other parameters in the basic reproduction number ${\mathcal{R}}_0$ be the same as those from Italy in Table 3.

Figure 5a gives the intersection line of the surface ${\mathcal{R}}_0={\mathcal{R}}_0({\mu }_3,{\gamma }_3)$ and the plane ${\mathcal{R}}_0=1$. Figure 5b shows the contour lines corresponding to the different values of ${\mathcal{R}}_0$. The contour lines give the range of the parameters, ${\mu }_3$ and ${\gamma }_3$, which are directly related to environmental and re-positive factors. Particularly, the area on the right side of the contour line with ${\mathcal{R}}_0=1$ can be seen as a high-risk area, which means that the disease infection will exist forever, and the area on the left side of the contour line with ${\mathcal{R}}_0=1$ can be seen as a low-risk area, which means that the disease infection will eventually be cleared. That is, when the parameter related to environmental factors (${\gamma }_3$) and the parameter related to re-positives (${\mu }_3$) have their values in the low-risk area, the disease infection is controllable.

Moreover, Figure 6a,b shows that when the parameters ${\gamma }_3$ (the infection rate due to environmental vectors) or ${\mu }_3$ (the rate of re-positives) decreases, the cumulative number of infected cases also decreases significantly.

According to the above figures, the re-positive rate ${\mu }_3$ decreased 10 times, and the cumulative number of infected cases decreased $48\%$. This is alarming, and suggests we need to pay more attention to recovered patients and their potential infectivity. We may need to re-evaluate hospitals’ discharge criteria and current patient management systems. Therefore, we stress that caution should be exercised even after recovery from SARS-CoV-2. We believe that all discharged patients should undergo medical observation and quarantine for at least 14 days. Longer periods of observation and surveillance may be necessary. Patients who have been infected with COVID-19 virus in the past should also comply with epidemiological control measures, such as wearing a mask and keeping distance from others.

The environment-related parameter ${\gamma }_3$ decreased 10 times, and the cumulative number of infected cases decreased $58\%$. However, some recent research has shown that there are many ways of communication for the novel coronavirus, and communication by environmental vectors is an important link that cannot be ignored. When an infected person sneezes or coughs, droplets are released that can contaminate surfaces. If a susceptible person touches these contaminated surfaces and then touches their mouth, nose, or eyes, they can become infected. If we want to ensure that the virus does not spread on a large scale, it is very important to sanitize the environment. For normal households, sanitizing with alcohol and disinfectant water can meet the demand, and families with susceptible conditions can also purchase ozone disinfection equipment for sanitizing. For places with local epidemic break-out, professional institutions are needed for sanitization. We have obtained some crucial epidemiological parameters that need greater emphasis for the mitigation and control of the COVID-19 pandemic.

4. Conclusions

Inspired by the diversity and complexity of COVID-19 infections, two key factors were considered in model (1). The first is that environmental vectors polluted by virus particles may lead to disease infection, and the second is that recovered individuals who have lost sufficient immune protection may return to the exposed class (i.e., called re-positive). In the past two years especially, these factors have been increasingly observed in the spread of COVID-19 infection, and are bringing great challenges and difficulties to the control of disease infection.

Theorems 1 and 2 give a complete characterization of the global dynamics of model (1). In other words, when the basic reproduction number ${\mathcal{R}}_0<1$, the disease-free equilibrium $Q_0$ is globally asymptotically stable, meaning that the disease infection will eventually be cleared. When the basic reproduction number ${\mathcal{R}}_0>1$, model (1) is uniformly persistent, meaning that the disease infection will exist forever and become an endemic disease. Furthermore, Theorem 3 gives some sufficient conditions for the global asymptotic stability of the endemic equilibrium $Q^{\ast }$ of model (1). It should be mentioned here that the sufficient conditions in Theorem 3 are very conservative and can be further improved.

Model (1) was applied to the data fitting of the cumulative confirmed cases and recovered cases of COVID-19 in Japan (20 May–18 June 2022) and Italy (12 January–10 February 2022), and good predictability of the data was observed.

Furthermore, based on data from Italy (12 January–10 February 2022), it was observed that the basic reproduction number ${\mathcal{R}}_0$ also had stronger dependence on the parameters ${\mu }_3$, ${\gamma }_2$, ${\gamma }_3$, and $\rho$. This indicates the necessity and practical significance of considering environmental and re-positive factors in the prevention and control of COVID-19 infection.