Analysis of a Model for Coronavirus Spread

Abstract

The spread of epidemics has always threatened humanity. In the present circumstance of the Coronavirus pandemic, a mathematical model is considered. It is formulated via a compartmental dynamical system. Its equilibria are investigated for local stability. Global stability is established for the disease-free point. The allowed steady states are an unlikely symptomatic-infected-free point, which must still be considered endemic due to the presence of asymptomatic individuals; and the disease-free and the full endemic equilibria. A transcritical bifurcation is shown to exist among them, preventing bistability. The disease basic reproduction number is calculated. Simulations show that contact restrictive measures are able to delay the epidemic’s outbreak, if taken at a very early stage. However, if lifted too early, they could become ineffective. In particular, an intermittent lock-down policy could be implemented, with the advantage of spreading the epidemics over a longer timespan, thereby reducing the sudden burden on hospitals.

1. Introduction

The coronavirus infection has been spreading for a few months. Authorities in several countries have relied on scientific tools for fighting the epidemics. With the lack of a vaccine, distancing methods have been forced on populations to avoid the transmission by direct contact. In laboratories, possible vaccines are being developed, but at the moment they are still at the experimental stage [1]. Meanwhile several models, mathematical, statistical and computer-science-based, are being developed to study the disease and contribute to fighting it.

Models for the spread of epidemics are classic, and an excellent presentation is [2]. In general, the total population is partitioned into at least two classes, susceptibles and infectives, with migrations from the former to the latter by disease propagation through direct or indirect contact, if the disease is transmissible. Additionally, if it can be overcome but causes relapses, the infected can become susceptible again, after maybe going through an intermediate class of being recovered. More sophisticated versions include quarantined and exposed individuals. Some of these classes will be considered also in the present study and illustrated in detail before the model formulation process.

In [3] the disease evolution forecast in several of the most affected countries is attempted, using for that purpose, parameter estimation techniques to calibrate the model. The involved compartments are susceptibles, asymptomatic individuals and symptomatic ones, which in turn are partitioned into reported and unreported cases. In [4] a simple SIRI model is considered, in which the recovered could still contribute to the disease spreading. The model is then extended to account for a possible vaccine, which, unfortunately, at present is not yet available, although several laboratories worldwide are trying to develop and test it, as mentioned above. In [5] a dynamic model for the diffusion of Covid-19 has been proposed. The transmission network is made by the bats–hosts–reservoir–people compartments; compare also [1]. As it amounts to about 14 differential equations and 25 parameters, it is rather complex. From it, the authors have obtained a simplified version, consisting of six compartments and 13 parameters. Then, the disease basic reproduction number has been calculated.

Our aim here is the mathematical analysis of a slightly modified version of the simpler model in [5]. The most important change accounts for the fact that asymptomatic people may indeed turn into fully symptomatic and infectious individuals. This feature also distinguishes the system introduced here from the one studied in [6], which, however, contains more compartments. The main aim of that study is the forecast of the epidemic spread in various cities in China, considering, additionally, weather data, which finally indicate that higher humidity favors the containment of a coronavirus epidemic. Our focus in the first part of this investigation is the theoretical analysis of the proposed system, and then we perform some preliminary simulations with realistic parameter values. More extended simulations will be devoted to a subsequent study.

The analysis of dynamical systems usually considers the possible equilibria that can be attained, assessing their feasibility and stability, and possible connections between them. For more details on these issues we refer the reader to classical texts, such as [7,8,9].

The paper is organized as follows. The main findings are outlined in the next section, which also discusses the results of numerical simulations. Section 3 contains an evaluation of their implications under various distancing policies. We formulate the model in Section 4, where we also analyze it mathematically, showing boundedness of the trajectories, establishing an expression for the disease basic reproduction number, finding its equilibria and assessing their local stability, and global stability is established just for the disease-free equilibrium. The section ends with the details on the numerical simulations.

2. Results

2.1. Theoretical Findings

The main analytical findings of this investigation are summarized in

Table 1. Equilibria of the model (1) and their feasibility conditions.

Equilibrium	Populations	Feasibility
$C_0$$=(S_0,0,0,0,0,0)$	$S_0={\displaystyle \frac{\Lambda }{d_p}}$	–
$C_I$$=(S_I,E_I,0,A_I,R_I)$	$S_I={\displaystyle \frac{\Lambda -B_T{E}_I}{d_p}}$
	$E_I={\displaystyle \frac{1}{B_T}}\left(\Lambda -{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I{\omega }_p^{\prime }{D}_T}}\right)$	$\Lambda >{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I{\omega }_p^{\prime }{D}_T}}$
	$A_I=\left({\displaystyle \frac{{\omega }_p^{\prime }}{H_T}}\right)E_I$	(resp. $R_0>1$)
	$R_I={\displaystyle \frac{{\gamma }_p^{\prime }}{d_p}}\left({\displaystyle \frac{{\omega }_p^{\prime }}{H_T}}\right)E_I$	$\alpha =1$ and $\xi =0$
$C^{\ast }$$=(S^{\ast },E^{\ast },I^{\ast },A^{\ast },R^{\ast })$	$S^{\ast }={\displaystyle \frac{\Lambda -B_T{E}^{\ast }}{d_p}}$	$\Lambda >$${\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}$
	$E^{\ast }={\displaystyle \frac{1}{B_T}}\left(\Lambda -{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right)$	(resp. $R_0>1$)
	$I^{\ast }=\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }\xi }{C_T{H}_T}}\right)E^{\ast }$	$\alpha \ne 1$ or $\xi \ne 0$
	$A^{\ast }=\left({\displaystyle \frac{\alpha {\omega }_p^{\prime }}{H_T}}\right)E^{\ast }$
	$R^{\ast }=\left[{\displaystyle \frac{{\gamma }_p}{d_p}}\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }\xi }{C_T{H}_T}}\right)\right.$$\left.+{\displaystyle \frac{{\gamma }_p^{\prime }}{d_p}}\left({\displaystyle \frac{\alpha {\omega }_p^{\prime }}{H_T}}\right)\right]E^{\ast }$

and

Table 2. Stability conditions of the equilibria of the model (1).

Point	Coefficients	Stability
$C_0$	$a_2=B_T+C_T+H_T$
	$a_1=H_T[B_T+C_T]+B_T{C}_T-{\beta }_I{S}_0((1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime })$	$\Lambda <{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}$
	$a_0=B_T{C}_T{H}_T-{\beta }_I{S}_0\left((1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right)$	(resp. $R_0<1$)
$C_I$	$b_3={\beta }_{I}k{A}_I+d_p+H_T+B_T+C_T$,
	$b_2=[{\beta }_{I}k{A}_I+d_p](H_T+B_T+C_T)$	$\Lambda >{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I{\omega }_p^{\prime }{D}_T}}$
	$+B_T{C}_T+H_T(B_T+C_T)-{\omega }_p^{\prime }k{\beta }_I{S}_I$
	$b_1=[{\beta }_{I}k{A}_I+d_p][B_T{C}_T+H_T(B_T+C_T)]$	(resp. $R_0>1$)
	$+H_T{B}_T{C}_T-{\omega }_p^{\prime }(k{d}_p+D_T){\beta }_I{S}_I$
	$b_0=[{\beta }_{I}k{A}_I+d_p]H_T{B}_T{C}_T-d_p{\omega }_p^{\prime }{D}_T{\beta }_I{S}^{\ast }$	$\alpha =1$ and $\xi =0$
$C^{\ast }$	$c_3={\beta }_I(I^{\ast }+k{A}^{\ast })+d_p+H_T+B_T+C_T$,
	$c_2=[{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p](H_T+B_T+C_T)$	$\Lambda >{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}$
	$+B_T{C}_T+H_T(B_T+C_T)-[\alpha {\omega }_p^{\prime }k+(1-\alpha ){\omega }_p]{\beta }_I{S}^{\ast }$	(resp. $R_0>1$)
	$c_1=[{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p][B_T{C}_T+H_T(B_T+C_T)]+H_T{B}_T{C}_T$	$\alpha \ne 1$ or $\xi \ne 0$
	$-[\alpha {\omega }_p^{\prime }(k{d}_p+D_T)+(1-\alpha ){\omega }_p(d_p+H_T)]{\beta }_I{S}^{\ast }$
	$c_0=[{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p]H_T{B}_T{C}_T$
	$-d_p[\alpha {\omega }_p^{\prime }{D}_T+(1-\alpha ){\omega }_p{H}_T]{\beta }_I{S}^{\ast }$

. The expressions of $B_T$, $C_T$, $H_T$, $D_T$ and $R_0$ are given by Equations (3) and (6).

The model (1) allows only three possible equilibria; the disease-free state $C_0$, where only susceptibles thrive; an equilibrium without symptomatic infected, which occurs only for a very particular case, when the exposed individuals become all asymptomatic infected; and finally, the endemic equilibrium $C^{\ast }$. All these equilibria are locally asymptotically stable, if suitable, rather complicated conditions, hold. Among the endemic and the disease-free equilibrium bistability is impossible, since they are related via a transcritical bifurcation.

2.2. Simulations Results

We have performed some simulations with the parameter values listed in

Table 3. Parameters values.

Parameter	Value	Parameter	Value
$\Lambda$	500	$d_p$	$8.2\times {10}^{-6}$
${\gamma }_p$	$1.764$	${\gamma }_p^{\prime }$	$0.6024$
$\xi$	$0.1$	k	$0.1\in$ [0:005; 0:2]
$\mu$	$0.001$	$\alpha$	$0.15\in [0.01,0.3]$
${\omega }_p$	$0.1$	${\omega }_p^{\prime }$	$0.1$
$\nu$	0

, to simulate various implementations of the distancing policy, which actually is in current use in several countries. The simulations may not be fully realistic, but our point is to investigate their qualitative behavior, not to give quantitative forecasts.

We look at the influence that the time of starting the restrictive measures has on the disease spread, while keeping fixed the time of their lifting. We next investigate the effect of the time at which the restricting measures are lifted.

Now comparing the results for the start of implementation at $t_1=1$ and $t_1=10$, and ending them at the same time, it is seen that the earlier the measures are taken, the better it is, because the epidemic’s outbreak is kept in check. In Figure 1 the epidemic outbreak starts around time 30, immediately after lifting the restrictions, while in Figure 2 the initial surge before the measures are implemented is damped by their implementation, and after their lifting the outbreak occurs. Both figures use $t_2=30$. The same result is seen using $t_2=90$ as the time for removing the restrictions; compare Figure 3 and Figure 4. In Figure 3 nothing apparently happens until time 100 because of the restrictions. When they are lifted, the epidemic spreads. In Figure 4 there is a small peak at the onset of the contagion, immediately curbed by the containment measures, lasting as long as they are in use. In spite of their longer implementation, the outbreak occurs nevertheless with the peak at the same time as in Figure 3.

The investigation of different timings for introducing and relaxing the distancing measures shows that a late implementation has no effect, as the peak of the epidemic occurs and then these measures are ineffective, independently of how long they are implemented. An earlier implementation followed by their subsequent lifting leads to a secondary peak at some time later, the occurrence of which seems to be related to the time for which the measures are implemented; the longer the latter, the longer the delay in the secondary outbreak. However, the number of affected people remains the same.

Unfortunately, the result of the simulations indicates that essentially the whole population gets affected by the disease. Only the timings differ, if distancing measures are taken.

Thus, if the measures are implemented too late, independently of the time at which they are removed, the outbreak occurs and their subsequent application becomes, therefore, irrelevant, as it cannot be kept in check any longer; compare Figure 2 and Figure 4. On the other hand, by implementing them at the early stages of the contagion process, the outbreak can be delayed, as long as these measures are implemented, as can be seen from Figure 1 and Figure 3. If they are lifted, the final results of the epidemic’s outbreak are essentially the same as if they were not at all implemented, in terms of the number of people being affected by the disease and with possible ultimate fatal consequences; compare the peaks of all the infected classes in Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5 also with the results in Figure 6 where no measures are taken to prevent the epidemic from spreading.

An intermittent lock-down policy, simulated as an alternative way of coping with the outbreak, might be important to render the burden on hospitalizations smaller, as it tends to spread the epidemic over a longer timespan.

For the particular situation in Italy, note that patient number 1 was diagnosed on 21 February, and the distancing measures in the area were in place starting the following days up to about two weeks later, and then extended to the whole country. Incidentally, patient number 0, the initial carrier of the disease, has never been identified, although there are some speculations. However, in the current news, it is reported that the virus was already circulating yet not known of in Northern Italy in January, which means that additional time had elapsed before the restrictions were applied.

Thus, apparently, these results are negative as for the possibility of containing the spread in the long run, in line with what is hinted in [10], with the exception of the intermittent distancing measures policy, which may spread the epidemic’s effects over longer timespans. However, there are some assumptions in the model that make it too crude, so that we plan a deeper subsequent study. In particular, here the results depend on homogeneous mixing, which for a large country is hardly the case. Secondly, this is a continuous model, for which the compartments are depleted only asymptotically. Thus it is not possible to prevent the class of infectives from vanishing in finite time, so that even a small residual in it would start the epidemic’s outbreak again. Therefore the somewhat negative results obtained might hopefully be better off in practice. Suitable modifications of the model along these lines will be the subject of a further investigation.

3. Discussion

We have investigated a simple model for the coronavirus pandemic. The steady states, apart from a symptomatic-infected-free point, which is unlikely to exist, are the disease-free equilibrium and the endemic state. The model differs from other current models that are being studied for a few features. From the simplified model that appears in [5], because our formulation contains less equations, it does not consider the viruses compartment, and above all, we allow disease-related mortality, which apparently is missing in the cited paper. Furthermore, we allow the progression of asymptomatic individuals to the class of fully symptomatic. This feature certainly distinguishes it also from [6], where asymptomatics recover or become diagnosed with the disease, but do not spread it any longer. In the present situation in Italy our assumption is very realistic.

There is no possibility of bistability in our situation, as the two fully meaningful equilibria are related to each other via a transcritical bifurcation. The disease-free equilibrium is also globally asymptotically stable, if it is locally asymptotically stable. An expression for the basic reproduction number is established, with a possibly realistic numerical value [11,12].

The simulations show that containment measures could be effective in delaying the epidemic’s outbreaks if taken at a very early stage, but when lifted the outbreaks would occur anyway and affect almost the whole population. However, this last statement should be mitigated by the drawbacks inherent in the model’s assumptions, as mentioned in the previous section, thereby leaving hope that in practice it will not occur, if the measures are properly implemented.

We next discuss in detail the various different restriction policies that we have simulated.

3.1. Epidemic with a Lock-Down

In this case, in particular, assuming for the disease transmission coefficient the reference value ${\beta }_I={10}^{-7}$, we reduce it to ${\beta }_I={10}^{-10}$ during the interval $[t_1,t_1+t_2]$ and reinstate the standard value afterwards; we monitor the epidemic’s evolution over six months. Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5 show the results of different choices for $t_1$ and $t_2$. Containment measures are effective as long as they are implemented, if they are taken early enough, before the epidemic attains its peak.

Since reducing the transmission by one order of magnitude means that to infect a susceptible with rate ${\beta }_I$, it is necessary for only one infected; with ${\beta }_I/10$, 10 infected would be necessary. Thus since the lock-down is not perfect, as for instance, some essential activities like food production are still going on, a hypothetical reasonable estimate for the contact reduction is three orders of magnitude. A comparison with a different, milder reduction, ${\beta }_I={10}^{-8}$ is made, showing essentially no difference in the results, see Figure 7.

3.2. Epidemic with Total Isolation

We changed also the policy to an improbable absolute confinement of every individual in the population, reducing the transmission to exactly zero. The results show no change with respect to those of the lock-down policy. We report only Figure 5, which is identical to Figure 1. The same occurs in the cases contemplated by Figure 2, Figure 3 and Figure 4.

3.3. The Simplified No-Demographics Model

We repeated the simulations for the model (1) in which we set $\Lambda =d_p=0$. In the simulations we observed some small changes in the susceptibles behavior, with respect to the full model with vital dynamics. Figure 8 and Figure 9 are the counterparts of the Figure 1 and Figure 2. The ultimate impact of the epidemic is essentially the same; compare in particular, the curves of recovered and deceased. For the total isolation case, Figure 10 shows the same features; compare it with Figure 5. Similar considerations hold for the various remaining cases, and therefore, the pictures are not reported.

3.4. Investigation of Different Timings for Restrictions’ Introduction and Lifting

A further study has been carried out to assess the impact of the time until taking action on the containment measures. All the possible different combinations of simple restriction or total isolation as well as the presence or the absence of demographic effects give essentially the same results. Therefore we present only the results for some selected alternatives, giving the plots in semilogarithmic or total population values, but stressing that for the options not considered, the figures would be the same.

In case the first restriction measure is taken too late, specifically at time $t=120$, and followed by lifting it either one month or three months later, the epidemic occurs and the measures have no effect whatsoever; see Figure 11, where measures are kept for three months.

Beginning the restrictions after three months from the start of the epidemic and removing them one month afterwards, causes a second peak about two months later; i.e., six months after the onset of the disease spreading (Figure 12), with a higher number of affected individuals. If instead the lock-down is implemented for three months, the second peak is delayed further, occurring about three months later, Figure 13. Although the pictures are shown on different population scales, absolute values and semilogarithmic, a comparison of the heights of the peaks for the various types of infected subpopulations indicates no difference. Hence, these policies cannot significantly influence the number of people ultimately affected by the disease.

3.5. The Intermittent Lock-Down Policy

We finally simulated a policy that attempts to assess the number of infectives at regular times, with a of period one week. If they exceed a threshold, taken to be 10, the lock-down is implemented for a week, and then lifted. Figure 14 and Figure 15 show the results for the case with vital dynamics and in the case of $\Lambda =d_p=0$. Note that susceptibles in both cases are at a constant value, the vertical scale being extremely small. The infected are kept below the threshold, and the periodic recurrences of the epidemic somewhat change its final impact, as the curves of recovered are reduced by about two orders of magnitude, and above all, the ones of the deceased decrease by about four orders, with respect to the ones found with the one-time lock-down policy. The other relevant change is that here the phenomenon is observed over a longer timespan. Thus the cumulative effects are spread out over a much longer time. This will have some importance to lessening the burden on hospitals. Figure 16 shows the results if the check policy starts immediately at time 1 rather than after a week.

Comparing the population values with the intermittent policy with the one time lock-down, done early enough and implemented for one month, the final outcomes are milder than the latter. Thus the intermittency allows the control of the outbreaks. Susceptibles are almost depleted in the one-time policy; with the intermittent one, however, they are essentially spared from getting the disease; compare Figure 17 and Figure 18.

The intermittency has also been checked with different time intervals. Comparing Figure 19, Figure 20, Figure 21 and Figure 22, it is seen that the more frequent the checks are implemented, the lower are the peaks in the exposed class, which in turn leads to a smaller cumulative number of recovered and fatalities, at least comparing the policies for the one- and two-weeks alternatives, Figure 19 and Figure 20. For the longer intervals between the checks, again the peaks are higher, the longer the timespan, but it is observed that as time elapses, their heights tend to decrease; see Figure 21 and Figure 22.

4. Materials and Methods

Here we develop a mathematical model of coronavirus, which is a zoonotic disease. Its biological characteristics indicate that the virus transmission occurred first from infected animals to humans [5], and then spread among populations worldwide by contact with infected individuals, to make it a pandemic.

Let $N\left(t\right)$ denote the total population. It is partitioned into the following five disjoint classes of individuals:

$S\left(t\right)$: The susceptible class, the individuals who have not yet been exposed to the virus.

$E\left(t\right)$: The exposed class, describing people who have become in contact with the virus, but are in incubation period and not yet able to spread the disease; possible presymptomatic individuals that can transmit the infection [13,14,15] are assumed to have already moved to the asymptomatic class defined below.

$I\left(t\right)$: The symptomatic infectious class, individuals that manifest symptoms and can spread the disease.

$A\left(t\right)$: The asymptomatic infectious class; those persons that can spread the disease even without having explicit symptoms.

$R\left(t\right)$: The removed class, that includes the people that recovered from the disease.

Thus, $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+R\left(t\right)$.

The basic mechanisms underlying the model are shown in Figure 23. The model is formulated taking into account all the possible interactions among the compartments that were described above.

Under the quasi-steady-state assumption of the total human population, we impose that susceptible individuals are recruited at the constant rate $\Lambda$, become infected by direct contact with an infectious individual at rate ${\beta }_I$, which is scaled by a factor k to account for the possibility that the latter is asymptomatic. Finally, all human individuals are subject to natural mortality $d_p$. These considerations are incorporated in the first equation of the system (1).

Individuals that contract the disease are accounted for in the second equation of (1). They become exposed, i.e., they cannot yet spread the virus, which needs an incubation period within the body of its hosts. In this class enter the susceptibles that were contaminated in the two ways described earlier. People leave it by becoming infectious, and either showing symptoms, thereby migrating into class I, or not, therefore, finding themselves in class A. The progression rates into these two classes are ${\omega }_p$ and ${\omega }_p^{\prime }$. Furthermore, we assume that a fraction $\alpha$ becomes asymptomatic and $1-\alpha$ instead will manifest symptoms.

The third equation models the symptomatic infectious, recruited from the exposed class at rate $(1-\alpha )\omega$ as described above. Furthermore, there could be asymptomatic individuals that become symptomatic at rate $\xi$. They leave this class by either progressing to the recovered class at rate ${\gamma }_p$, or dying, naturally or by causes related to the disease at rate $\mu$.

The asymptomatic individuals modeled in the fourth equation appear from the exposed ones, and leave the class by overcoming the disease at rate ${\gamma }_p^{\prime }$, dying naturally or by disease-related causes at rate $\nu$, or eventually showing the symptoms, for which they migrate into class I.

Recovered individuals are those that have healed from the disease. They are subject only to natural mortality. We assume that they have also become immune so that they are unaffected if become in contact with the infectious.

Note that in the simulations also the cumulative class of disease-related deceased people is shown, although the dead are not explicitly accounted for in the model. They indeed represent a sink, and thus do not contribute to the disease propagation. Incidentally, instead, in cultures where the deceased are kept for a while before burial and become in contact with the relatives, it may be necessary to introduce this class in the model, as another potential source of infection.

Taking into account the above considerations, the model dynamics is regulated by the following system of nonlinear ordinary differential equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \Lambda -{\beta }_{I}S(I+kA)-d_{p}S,\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& =& {\beta }_{I}S(I+kA)-(1-\alpha ){\omega }_{p}E-\alpha {\omega }_p^{\prime }E-d_{p}E,\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =& (1-\alpha ){\omega }_{p}E-({\gamma }_p+d_p+\mu )I+\xi A,\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =& \alpha {\omega }_p^{\prime }E-({\gamma }_p^{\prime }+d_p+\nu )A-\xi A,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =& {\gamma }_{p}I+{\gamma }_p^{\prime }A-d_{p}R,\hfill \end{array}$$

(1)

or alternatively, excluding completely the demographic features, by setting $\Lambda =0$ and $d_p=0$ in (1). All the parameters are nonnegative and their meaning is summarized in

Table 4. Model parameters and their meaning.

$\Lambda$	susceptibles recruitment rate
$d_p$	natural mortality
${\beta }_I$	disease transmission rate
k	transmissibility ratio between asymptomatics and symptomatics
$\mu$	disease-related mortality for infected
$\nu$	disease-related mortality for asymptomatics
${\omega }_p$	progression rate from exposed to symptomatic
${\omega }_p^{\prime }$	progression rate from exposed to asymptomatic
$\alpha$	fraction of exposed that turn asymptomatic
$\xi$	progression rate from asymptomatic to symptomatic
${\gamma }_p$	recovery rate from symptomatic infection
${\gamma }_p^{\prime }$	recovery rate from asymptomatic infection

. Note that in view of the definitions,

$$\frac{1}{{\omega }_p},\frac{1}{{\omega }_p^{\prime }}\frac{1}{{\gamma }_p},\frac{1}{{\gamma }_p^{\prime }},$$

represent respectively the incubation period before manifesting symptoms, the latent period before becoming asymptomatic infectious, the infectious period for symptomatic infection and the infectious period for asymptomatic infection.

Theorem 1.

The system trajectories are bounded. Letting

$$M=max\left\{N\left(0\right),{\displaystyle \frac{\Lambda }{d_p}}\right\}$$

the set

$$\begin{array}{c}\hfill \Gamma =\left\{\right(S,E,I,A,R):S+E+I+A+R\le M,S>0,E\ge 0,I\ge 0,A\ge 0,R\ge 0\}.\end{array}$$

(2)

represents their ultimate attractor. In particular, if $N\left(0\right)<\Lambda d_p^{-1}$, $M=\Lambda d_p^{-1}$.

Proof. 

From the system (1) it follows that the total population evolves as follows:

$${\displaystyle \frac{dN}{dt}}+d_{p}N=\Lambda -\nu A-\mu I\le \Lambda .$$

Solving the differential inequality easily gives

$$N\left(t\right)\le N\left(0\right)exp(-d_{p}t)+{\displaystyle \frac{\Lambda }{d_p}}[1-exp(-d_{p}t)]\le M,$$

so that all subpopulations, being nonnegative, are bounded as well. □

Note that $\Gamma$ is positively invariant since all solutions of system (1) originating in $\Gamma$ remain there for all $t>0$, in view of the existence and uniqueness of its solutions.

4.1. System’s Equilibria Assessment

The equilibrium points of the model are obtained by equating the right hand side of system (1) to zero. The solution of the so-obtained algebraic system gives three equilibrium points: the coronavirus-free equilibrium $C_0=(S_0,0,0,0,0,)$, the coronavirus-symptomatic-infected-free equilibrium $C_I=(S_I,E_I,0,A_I,R_I)$ with conditions $\alpha =1$ and $\xi =0$, and the fully coronavirus endemic equilibrium $C^{\ast }=(S^{\ast },E^{\ast },I^{\ast },A^{\ast },R^{\ast })$ when either $\alpha \ne 1$ or $\xi \ne 0$. Specifically, for the former two we have:

$$S_0={\displaystyle \frac{\Lambda }{d_p}},E_I={\displaystyle \frac{1}{B_T}}\left(\Lambda -{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I{\omega }_p^{\prime }{D}_T}}\right),S_I={\displaystyle \frac{\Lambda -B_T{E}^{\ast }}{d_p}},A_I=\left({\displaystyle \frac{{\omega }_p^{\prime }}{H_T}}\right)E_I,R_I=\left({\displaystyle \frac{{\gamma }_p^{\prime }}{d_p}}{\displaystyle \frac{{\omega }_p^{\prime }}{H_T}}\right)E_I,$$

where

$$B_T=(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }+d_p,C_T={\gamma }_p+\mu +d_p,D_T=\xi +k({\gamma }_p+\mu +d_p),H_T={\gamma }_p^{\prime }+\nu +\xi +d_p.$$

(3)

The feasibility conditions for $C_I$ are

$$\Lambda >{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I{\omega }_p^{\prime }{D}_T}},\alpha =1\mathrm{and}\xi =0.$$

(4)

For the fully endemic equilibrium we find

$$\begin{array}{cc}{E}^{\ast }={\displaystyle \frac{1}{B_T}}\left(\Lambda -{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right),\hfill & S^{\ast }={\displaystyle \frac{\Lambda -B_T{E}^{\ast }}{d_p}},\hfill \\ I^{\ast }=\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }\xi }{C_T{H}_T}}\right)E^{\ast },\hfill & A^{\ast }=\left({\displaystyle \frac{\alpha {\omega }_p^{\prime }}{H_T}}\right)E^{\ast }\hfill \\ \mathrm{and} R^{\ast }=\left({\displaystyle \frac{{\gamma }_p}{d_p}}\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }\xi }{C_T{H}_T}}\right)+{\displaystyle \frac{{\gamma }_p^{\prime }}{d_p}}\left({\displaystyle \frac{\alpha {\omega }_p^{\prime }}{H_T}}\right)\right)E^{\ast },\hfill \end{array}$$

with feasibility condition

$$\Lambda >{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}},\mathrm{and} \mathrm{either}\alpha \ne 1or\xi \ne 0.$$

(5)

4.2. The Basic Reproduction Number

The basic reproduction number $R_0$ for system (1) is found using the next generation matrix method [16]. The reduced system of (1) may be written in compact form as: $X^{\prime }=F\left(X\right)-V\left(X\right)$ where $X=(E,I,A)$

$$F(E,I,A)=\left(\begin{array}{c}{\beta }_{I}S(I+kA)\hfill \\ 0\hfill \\ 0\hfill \end{array}\right),V(E,I,A)=\left(\begin{array}{c}-(1-\alpha ){\omega }_{p}E-\alpha {\omega }_p^{\prime }E-d_{p}E\hfill \\ (1-\alpha ){\omega }_{p}E-({\gamma }_p+d_p+\mu )I+\xi A\hfill \\ \alpha {\omega }_p^{\prime }E-({\gamma }_p^{\prime }+d_p+\nu )A-\xi A\hfill \end{array}\right).$$

The Jacobian matrices of $F\left(X\right)$ and $V\left(X\right)$ at the disease-free equilibrium point $C_0$ are

$$J_F\left(C_0\right)=\left(\begin{array}{ccc}0& {\beta }_I{S}_0& {\beta }_I{S}_{0}k\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$$

and

$$J_V\left(C_0\right)=\left(\begin{array}{ccc}-B_T& 0& 0\\ (1-\alpha ){\omega }_p& -C_T& \xi \\ \alpha {\omega }_p^{\prime }& 0& -H_T\end{array}\right).$$

We find that

$$J_V^{-1}\left(C_0\right)=\left(\begin{array}{ccc}{\displaystyle \frac{-1}{B_T}}& 0& 0\\ {\displaystyle \frac{-[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }\xi ]}{C_T{B}_T{H}_T}}& {\displaystyle \frac{-1}{C_T}}& {\displaystyle \frac{-\xi }{C_T{H}_T}}\\ {\displaystyle \frac{-\alpha {\omega }_p^{\prime }}{B_T{H}_T}}& 0& {\displaystyle \frac{-1}{H_T}}\end{array}\right).$$

The next generation matrix is

$$-J_F\left(C_0\right)J_V^{-1}\left(C_0\right)=\left(\begin{array}{ccc}{\displaystyle {\beta }_I{S}_0\frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{B}_T{H}_T}}& {\displaystyle \frac{{\beta }_I{S}_0}{C_T}}& {\displaystyle \frac{{\beta }_I{S}_0{D}_T}{C_T{H}_T}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

Thus

$$R_0=\rho (-J_F\left(C_0\right)J_V^{-1}\left(C_0\right))={\beta }_I{S}_0\frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{B}_T{H}_T}.$$

(6)

The conditions (4) (resp. (5)) are equivalent to $R_0>1$ for $\alpha =1$ and $\xi =0$ (resp.$R_0>1$ for either $\alpha \ne 1$ or $\xi \ne 0$).

We have the following theorem

Theorem 2.

System (1) has the following equilibria:

1. 

The coronavirus-free equilibrium $C_0=(S_0,0,0,0,0)=\left(\frac{\Lambda }{d_p},0,0,0,0\right)$ which exists always.

2. 

In addition, if $R_0>0$ then system (1) admits another nontrivial equilibrium, in fact:

When $\alpha =1$ and $\xi =0$, it is the coronavirus-symptomatic-infected-free equilibrium $C_I=(S_I,E_I,I_I,A_I,R_I)$.

When either $\alpha \ne 1$ or $\xi \ne 0$, it is the fully coronavirus endemic equilibrium $C^{\ast }=(S^{\ast },E^{\ast },I^{\ast },A^{\ast },R^{\ast })$.

4.3. System’s Equilibria Stability

4.3.1. Local Stability

In this subsection we investigate the local stability of the system’s equilibria.

Theorem 3.

Letting

$$\begin{array}{cc}\hfill & a_2=B_T+C_T+H_T,\hfill \\ \hfill & a_1=H_T[B_T+C_T]+B_T{C}_T-{\beta }_I{S}_0((1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime })\hfill \\ \hfill & a_0=B_T{C}_T{H}_T-{\beta }_I{S}_0\left((1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right).\hfill \end{array}$$

(7)

1. 

The coronavirus-free equilibrium $C_0=(S_0,0,0,0,0)$ of the system (1) is locally asymptotically stable if

$$\Lambda <{\displaystyle \frac{d_p}{{\beta }_I}}{\displaystyle \frac{B_T{C}_T{H}_T}{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}},(\mathrm{resp}.R_0<1).$$

(8)

2. 

If $\Lambda >{\displaystyle \frac{d_p}{{\beta }_I}}{\displaystyle \frac{B_T{C}_T{H}_T}{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}}$, (resp. $R_0>1$), then the coronavirus-free equilibrium $C_0$ of the system (1) is unstable.

Proof. 

The Jacobian matrix of system (1) at the coronavirus-free equilibrium $C_0$ is:

$$J\left(C_0\right)=\left(\begin{array}{ccccc}-d_p& 0& -{\displaystyle \frac{{\beta }_I\Lambda }{d_p}}& -{\displaystyle \frac{k{\beta }_I\Lambda }{d_p}}& 0\\ 0& -B_T& {\displaystyle \frac{{\beta }_I\Lambda }{d_p}}& {\displaystyle \frac{k{\beta }_I\Lambda }{d_p}}& 0\\ 0& (1-\alpha ){\omega }_p& -C_T& \xi & 0\\ 0& \alpha {\omega }_p^{\prime }& 0& -H_T& 0\\ 0& 0& {\gamma }_p& {\gamma }_p^{\prime }& -d_p\end{array}\right).$$

At point $C_0$, the eigenvalues of J are $-d_p$ of multiplicity order two and the roots of the following characteristic polynomial of a three by three submatrix of J whose coefficients $a_i$, $i=0,\dots ,2$ are given in (7):

$${\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0.$$

(9)

It is evident that $a_2>0$. From condition (8) the following inequalities are also satisfied

$$\begin{array}{ccc}\hfill a_0& =& B_T{C}_T{H}_T-{\beta }_I{S}_0\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\hfill \\ & =& \frac{{\beta }_I}{d_p}\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\left(\frac{d_p}{{\beta }_I}\frac{B_T{C}_T{H}_T}{\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}-\Lambda \right)>0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]a_1& =& [H_T(B_T+C_T)+B_T{C}_T]\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\hfill \\ & & -{\beta }_I{S}_0\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right][(1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime }]\hfill \\ & =& [H_T(B_T+C_T)+B_T{C}_T]\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\hfill \\ & & +[a_0-B_T{C}_T{H}_T][(1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime }]\hfill \\ & =& H_T(B_T+C_T)(1-\alpha ){\omega }_p{H}_T+(H_T{C}_T+B_T{C}_T)\alpha {\omega }_p^{\prime }\xi \hfill \\ & & +a_0[(1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime }]>0\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]a_1{a}_2& =& H_T(B_T+C_T)(1-\alpha ){\omega }_p{H}_T{a}_2+(H_T{C}_T+B_T{C}_T)\alpha {\omega }_p^{\prime }\xi a_2\hfill \\ & & +a_0[(1-\alpha ){\omega }_p+k\alpha {\omega }_p^{\prime }]a_2\hfill \\ & =& H_T(B_T+C_T)(1-\alpha ){\omega }_p{H}_T{a}_2+(H_T{C}_T+B_T{C}_T)\alpha {\omega }_p^{\prime }\xi a_2\hfill \\ & & +a_0(1-\alpha ){\omega }_p(B_T+C_T+H_T)+a_{0}k\alpha {\omega }_p^{\prime }(B_T+C_T+H_T)\hfill \\ & >& B_T{C}_T\alpha {\omega }_p^{\prime }\xi a_2+a_0(1-\alpha ){\omega }_p{H}_T+a_0\alpha {\omega }_p^{\prime }(k{C}_T+\xi )-a_0\alpha {\omega }_p^{\prime }\xi \hfill \\ & >& B_T{C}_T{H}_T\alpha {\omega }_p^{\prime }\xi +a_0[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T]-B_T{C}_T{H}_T\alpha {\omega }_p^{\prime }\xi \hfill \\ & =& a_0[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T].\hfill \end{array}$$

Thus, $a_i>0$, $i=0,\dots ,2$ and $a_2{a}_1>a_0$.

Then, according to the Routh–Hurwitz criterion, all the roots of the characteristic Equation (9) have negative real parts. Therefore, the coronavirus-free equilibrium point $C_0$ is locally asymptotically stable under condition (8). □

Since we can deduce the stability of the coronavirus symptomatic infected-free equilibrium $C_I$ from the stability of the coronavirus endemic equilibrium $C^{\ast }$ simply by taking $\alpha =1$ and $\xi =0$ in the latter, we now just analyze the coronavirus endemic equilibrium $C^{\ast }$.

Theorem 4.

Let

$$\left\{\begin{array}{ccc}{c}_3& =& {\beta }_I(I^{\ast }+k{A}^{\ast })+d_p+H_T+B_T+C_T>0,\hfill \\ c_2& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p](H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & -[\alpha {\omega }_p^{\prime }k+(1-\alpha ){\omega }_p]{\beta }_I{S}^{\ast },\hfill \\ c_1& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p][B_T{C}_T+H_T(B_T+C_T)]+H_T{B}_T{C}_T\hfill \\ & & -[\alpha {\omega }_p^{\prime }(k{d}_p+D_T)+(1-\alpha ){\omega }_p(d_p+H_T)]{\beta }_I{S}^{\ast },\hfill \\ c_0& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p]H_T{B}_T{C}_T-d_p[\alpha {\omega }_p^{\prime }{D}_T+(1-\alpha ){\omega }_p{H}_T]{\beta }_I{S}^{\ast }.\hfill \end{array}\right.$$

(10)

The coronavirus endemic equilibrium $C^{\ast }$ is locally asymptotically stable if

$$\Lambda >{\displaystyle \frac{d_p}{{\beta }_I}}{\displaystyle \frac{B_T{C}_T{H}_T}{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}},(resp.R_0>1).$$

(11)

Proof. 

The Jacobian matrix of system (1) at the coronavirus endemic equilibrium $C^{\ast }$ is:

$$J\left(C^{\ast }\right)=\left(\begin{array}{ccccc}-{\beta }_I(I^{\ast }+k{A}^{\ast })-d_p& 0& -{\beta }_I{S}^{\ast }& -{\beta }_I{S}^{\ast }k& 0\\ {\beta }_I(I^{\ast }+k{A}^{\ast })& -B_T& {\beta }_I{S}^{\ast }& {\beta }_I{S}^{\ast }k& 0\\ 0& (1-\alpha ){\omega }_p& -C_T& \xi & 0\\ 0& \alpha {\omega }_p^{\prime }& 0& -H_T& 0\\ 0& 0& {\gamma }_p& {\gamma }_p^{\prime }& -d_p\end{array}\right).$$

At point $C^{\ast }$, the eigenvalues of J are $-d_p$ and the roots of the characteristic polynomial of a three by three submatrix of J. The characteristic equation, in which the coefficients $c_i$, $i=0,\dots ,3$ are given in (10), is:

$${\lambda }^4+c_3{\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0.$$

(12)

It is evident that $c_3>0$. From condition (11) the following inequalities are also satisfied.

$$\begin{array}{ccc}\hfill c_0& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p]H_T{B}_T{C}_T-{\beta }_I[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T]d_p{S}^{\ast }\hfill \\ & =& {\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }(\xi +k{C}_T)\right]B_T{E}^{\ast }+d_p{H}_T{C}_T{B}_T\hfill \\ & & -{\beta }_I[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T](\Lambda -B_T{E}^{\ast })\hfill \\ & =& {\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\left(2B_T{E}^{\ast }+\frac{d_p{H}_T{C}_T{B}_T}{{\beta }_I[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T]}-\Lambda \right)\hfill \\ & =& {\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]B_T{E}^{\ast }>0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill c_1& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p][B_T{C}_T+H_T(B_T+C_T)]+H_T{B}_T{C}_T\hfill \\ & & -[(1-\alpha ){\omega }_p(d_p+H_T)+\alpha {\omega }_p^{\prime }(k{d}_p+D_T)]{\beta }_I{S}^{\ast }\hfill \\ & =& \left[{\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)E^{\ast }+d_p\right][B_T{C}_T+H_T(B_T+C_T)]\hfill \\ & & -{\beta }_I[(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }k](\Lambda -B_T{E}^{\ast })\hfill \\ & =& {\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T(B_T+C_T)}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }\xi B_T}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }{D}_T(B_T+H_T)}{H_T}}\right)E^{\ast }\hfill \\ & & +d_p[B_T{C}_T+H_T(B_T+C_T)]-{\beta }_I[(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }k](\Lambda -2B_T{E}^{\ast })\hfill \\ & =& {\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T(B_T+C_T)}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }\xi B_T}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }{D}_T(B_T+H_T)}{H_T}}\right)E^{\ast }\hfill \\ & & +d_p[B_T{C}_T+H_T(B_T+C_T)]-{\beta }_I[(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }k]\left({\displaystyle \frac{2d_p{B}_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}-\Lambda \right)\hfill \\ & =& {\beta }_I[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T]\left({\displaystyle \frac{C_T{H}_T+B_T(C_T+H_T)}{C_T{H}_T}}\right)E^{\ast }\hfill \\ & & +d_p\left({\displaystyle \frac{(B_T+C_T)H_T^2(1-\alpha ){\omega }_p+H_T{B}_T\alpha {\omega }_p^{\prime }\xi +C_T(B_T+H_T)\alpha {\omega }_p^{\prime }{D}_T}{\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right)>0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill c_2& =& [{\beta }_I(I^{\ast }+k{A}^{\ast })+d_p](H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & -[\alpha {\omega }_p^{\prime }k+(1-\alpha ){\omega }_p]{\beta }_I{S}^{\ast }\hfill \\ & =& \left[{\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)E^{\ast }+d_p\right](H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & -[\alpha {\omega }_p^{\prime }k+(1-\alpha ){\omega }_p]{\beta }_I{\displaystyle \frac{\Lambda -B_T{E}^{\ast }}{d_p}}\hfill \\ & =& {\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p(H_T+B_T)}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }\xi (H_T+B_T+C_T)}{C_T{H}_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }k(B_T+C_T)}{H_T}}\right)E^{\ast }\hfill \\ & & +d_p(H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & +{\beta }_I[(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }k]\left({\displaystyle \frac{(d_p+B_T)E^{\ast }-\Lambda }{d_p}}\right)\hfill \\ & =& {\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p(H_T+B_T)}{C_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }\xi (H_T+B_T+C_T)}{C_T{H}_T}}+{\displaystyle \frac{\alpha {\omega }_p^{\prime }k(B_T+C_T)}{H_T}}\right)E^{\ast }\hfill \\ & & +d_p(H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & +\frac{{\beta }_I[(1-\alpha ){\omega }_p+\alpha {\omega }_p^{\prime }k]}{B_T}\left(\Lambda -{\displaystyle \frac{(d_p+B_T)B_T{C}_T{H}_T}{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right)\hfill \\ & =& {\beta }_I(B_T+C_T+H_T)\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)E^{\ast }\hfill \\ & & +d_p(H_T+B_T+C_T)+C_T{H}_T+B_T\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T^2+H_T\alpha {\omega }_p^{\prime }\xi +\alpha {\omega }_p^{\prime }{C}_T{D}_T}{\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right)>0\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill c_1(c_3{c}_2-c_1)& =& {\beta }_I\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)E^{\ast }{c}_1{c}_2\hfill \\ & & +{\beta }_I\left[(H_T+C_T+d_p)H_T+C_T(C_T+d_p)\right]\left({\displaystyle \frac{\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}{C_T{H}_T}}\right)E^{\ast }{c}_1\hfill \\ & & +{\beta }_I(H_T+B_T+C_T+d_p)B_T\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)E^{\ast }{c}_1\hfill \\ & & +(H_T+B_T+C_T)(H_T+B_T+C_T+d_p)d_p{c}_1+C_T{H}_T(C_T+H_T+B_T)c_1\hfill \\ & & +B_T(B_T+C_T+H_T)\left({\displaystyle \frac{[(1-\alpha ){\omega }_p{H}_T^2+\alpha {\omega }_p^{\prime }{C}_T{D}_T]+\alpha {\omega }_p^{\prime }\xi H_T}{\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]}}\right)c_1\hfill \\ & >& {\beta }_I^2{B}_T{E}^{\ast 2}\left\{{\beta }_I\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]{\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)}^2{E}^{\ast }\right.\hfill \\ & & \left.+2(d_p+H_T+B_T+C_T)\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]\left({\displaystyle \frac{(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T}{C_T{H}_T}}\right)\right\}\hfill \\ & & +{\beta }_I{B}_T{(d_p+H_T+B_T+C_T)}^2\left[(1-\alpha ){\omega }_p{H}_T+\alpha {\omega }_p^{\prime }{D}_T\right]E^{\ast }\hfill \\ & =& c_0{c}_3^2.\hfill \end{array}$$

Thus, $c_i>0$, $i=0,\dots ,3$ and $c_1(c_3{c}_2-c_1)>c_0{c}_3^2$. Then, according to the Routh–Hurwitz criterion, all the roots of the characteristic Equation (12) have negative real parts. Therefore, the coronavirus endemic equilibrium point $C^{\ast }$ is locally asymptotically stable under condition (11). □

From Theorem 4 the following result is reached.

Theorem 5.

Let

$$\left\{\begin{array}{ccc}{b}_3& =& {\beta }_{I}k{A}_I+d_p+H_T+B_T+C_T>0,\hfill \\ b_2& =& [{\beta }_{I}k{A}_I+d_p](H_T+B_T+C_T)+B_T{C}_T+H_T(B_T+C_T)\hfill \\ & & -{\omega }_p^{\prime }k{\beta }_I{S}_I,\hfill \\ b_1& =& [{\beta }_{I}k{A}_I+d_p][B_T{C}_T+H_T(B_T+C_T)]+H_T{B}_T{C}_T\hfill \\ & & -{\omega }_p^{\prime }(k{d}_p+D_T){\beta }_I{S}_I,\hfill \\ b_0& =& [{\beta }_{I}k{A}_I+d_p]H_T{B}_T{C}_T-d_p{\omega }_p^{\prime }{D}_T{\beta }_I{S}^{\ast }.\hfill \end{array}\right.$$

(13)

The coronavirus symptomatic-infected-free equilibrium $C_I$ of the system (1) is locally asymptotically stable if

$$\Lambda >{\displaystyle \frac{d_p}{{\beta }_I}}{\displaystyle \frac{B_T{C}_T{H}_T}{{\omega }_p^{\prime }{D}_T}},(resp.R_0>1).$$

(14)

Proof. 

The result can easily obtained from Theorem 4 by taking $\alpha =1$ and $\xi =0$. □

Additionally, from the previous discussion, we can claim the following result:

Theorem 6.

There is a transcritical bifurcation between $C_0$ and $C^{\ast }$.

4.3.2. Global Stability

Next, we address the issue of global stability of the coronavirus–free equilibrium, employing as a tool a suitably constructed Lyapunov function and La Salle’s Invariance Principle.

Theorem 7.

The coronavirus-free equilibrium $C_0$ of model (1) is globally asymptotically stable if

$$\Lambda <{\displaystyle \frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}},(resp.R_0<1).$$

(15)

Proof. 

First, the four equations of (1) are independent of R, therefore, the last equation of (1) can be omitted without loss of generality. Hence, let us consider the following function:

$$P=\frac{1}{2S_0}{(S-S_0)}^2+E+\frac{B_T}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}(H_{T}I+D_{T}A)$$

(16)

It is easily seen that the above function is nonnegative and also $P=0$ if and only if $S=S_0$, $E=0$, $I=0$ and $A=0$. Further, calculating the time derivative of P along the positive solutions of (1), we find:

$$\begin{array}{ccc}\hfill \frac{dP}{dt}& =& \frac{1}{S_0}(S-S_0)(-{\beta }_{I}S(I+kA)-d_p(S-S_0))+{\beta }_{I}S(I+kA)-B_{T}E\hfill \\ & & +\frac{B_T{H}_T((1-\alpha ){\omega }_{p}E-C_{T}I+\xi A)}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}+\frac{B_T{D}_T(\alpha {\omega }_p^{\prime }E-H_{T}A)}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}\hfill \\ & =& -\frac{d_p}{S_0}{(S-S_0)}^2+{\beta }_I[2S-\frac{S^2}{S_0}-S_0](I+kA)+{\beta }_I{S}_0(I+kA)\hfill \\ & & -\frac{B_T{C}_T{H}_{T}I}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}-\frac{B_T(-\xi +D_T)H_{T}A}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}\hfill \\ & =& -\frac{d_p}{S_0}{(S-S_0)}^2+{\beta }_I[2S-\frac{S^2}{S_0}-S_0](I+kA)\hfill \\ & & +\left({\beta }_I{S}_0-\frac{B_T{C}_T{H}_T}{[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}\right)(I+kA)\hfill \\ & =& -\frac{d_p}{S_0}{(S-S_0)}^2+{\beta }_I[2S-\frac{S^2}{S_0}-S_0](I+kA)\hfill \\ & & +\frac{{\beta }_I}{d_p}\left(\Lambda -\frac{d_p{B}_T{C}_T{H}_T}{{\beta }_I[(1-\alpha ){\omega }_p{H}_T+D_T\alpha {\omega }_p^{\prime }]}\right)(I+kA).\hfill \end{array}$$

From condition (15) we can show that the coefficients of the term $I+kA$ in the last equality are negative. Further, we have $2S-\frac{S^2}{S_0}-S_0=-\frac{S^2-2S{S}_0+S_0^2}{S_0}=-\frac{{(S-S_0)}^2}{S_0}\le 0$ for all $S\ge 0$. Thus, we have $\frac{dP}{dt}\le 0$ for all $(S,E,I,A)\in {\mathbb{R}}_{+}^4$ and $\frac{dP}{dt}=0$ if and only if $(S,E,I,A)=(S_0,0,0,0)$. Thus, the only invariant set contained in ${\mathbb{R}}_{+}^4$ is $\left\{(S_0,0,0,0)\right\}$. Hence, La Salle’s theorem implies convergence of the solutions $(S,E,I,A)$ to $(S_0,0,0,0)$. From the last equation if (1) we can show obviously that R converge also to 0. Therefore $C_0$ is globally asymptotically stable if $R_0<1$. □

4.4. Numerical Simulations

The calculation of the value of $R_0$ according to (6) with the parameter values used in the numerical simulations gives $R_0=3.1402$, in line with the current estimates [11,12].

4.4.1. Simulations Methodology

We use a simple own-developed driver code calling the Matlab intrinsic routine ode45, implementing the classical Runge–Kutta 45 integration method for ordinary differential equations.

At first, we consider only the demographic simulation and show that the population is essentially at the same level during a year. This fact is substantiated also by the simulation results, for which there is scant difference between those of the model (1) and the ones obtained by using its no-demographic counterpart, where $\Lambda$ and $d_p$ are both set to zero.

We then perform three sets of simulations describing different possible scenarios. The first one considers lock-down, i.e., decreasing the contact rate significantly, but not to zero, as some essential activities are still open. Then the total isolation policy, for which the contact rate is set to zero. Finally an intermittent closure policy, for which when infectives reappear in a significant way, temporary lock-down measures resume again.

4.4.2. Data Acquisition

We use data published on official websites about the epidemic’s spread in Italy collected between 29 January and 28 March 2020, a period that spans 61 days, incremented by more recent information [17].

Using the day as the base time unit, we assume that the average incubation period lies in the interval between two and 14 days, with a mean of 8 days. Based on the percentage of the reported symptomatic infected patients, the proportion of symptomatic in the infected class $\alpha$ is estimated to be in the interval $[0.01,0.3]$. The correction k for asymptomatics to diffuse the disease is set in the range $k\in [0:005;0:2]$. There have been 27,359 deaths between 15 February and 29 April [17], with changes in the number of fatalities every day. Dividing the fatal cases by the timespan, one gets 370 daily fatalities, which gives a rate $0.0027$. Using this value in the simulation, puts the total losses to about ${10}^5$. But we observed that apparently children hardly get the disease, the younger and adult people have it generally in a mild form and fatalities occur mainly for the elderly people, compare with Figure 3 of [18]. In view of the fact that there is no age structure in this model, we corrected this value by taking a third of the above result to set the disease mortality rate at the final value $\mu =0.001$, which gives a reasonable estimate for the losses in the timespan, in rough agreement with the actual tallies. We neglect altogether mortality for the asymptomatics, setting $\nu =0$. Based on the officially published data we estimate ${\gamma }_p=0.1764$, ${\gamma }_p^{\prime }=0.6024$. For the initial values, the total population is obtained from the report published by the official cite of worldometers [19], $S=60461826$. To avoid demographic effects, we set the susceptible recruitment rate $\Lambda$ in order that on the timespan of the simulation the total population N does not change much.

4.4.3. The Pure Demographic Case

We simulate first the population model without disease. In so doing, we varied the parameter $\Lambda$ until a satisfactory behavior of N, the total population was found. With $\Lambda =500$ there is little variation of N during a whole year, the population remains roughly stable around the level $\mathrm{60,400,000}$, see Figure 24. In this way the demographic effects are sort of removed, and we can concentrate mainly on the epidemics. Actually, the number of newborns per day in Italy would be about four times higher, but as mentioned, we just would like here to hide the demographics from the simulations and not have a picture more adherent to reality.

4.4.4. Epidemics Spread in the Absence of Measures

Here we introduced the disease, with incidence ${\beta }_I={10}^{-7}$. The result is shown in Figure 25 for absolute numbers, and in Figure 6 in semilogarithmic scale. In this case no measures are assumed to be taken to counteract the epidemics. These results are reported in order be able to compare the simulations with restrictions to what would happen if the containment measures were not taken.

4.4.5. Containment Measures for the Epidemics

Finally, we considered the introduction of the distancing policy. It is assumed to start at time $t_1$ and end at time $t_1+t_2$. Two forms of containment measures are considered, substantially reducing the contact rate, or even setting it equal to zero, meaning the extreme measure of total individuals isolation.

In particular, we present the experience of using the reference value of the contact rate ${\beta }_I={10}^{-7}$, then reducing it to ${\beta }_I={10}^{-10}$ during the interval $[t_1,t_1+t_2]$. We then reset it to its previous reference value after time $t_1+t_2$. We monitored the epidemics evolution over six months.

The alternative, milder choice ${\beta }_I={10}^{-8}$ is also used, for comparison.

The simulations are then repeated with total isolation, setting ${\beta }_I=0$ during the implementation of the restrictions.

A comparison of the results with the model obtained by disregarding the demographic parameters, i.e., setting $\Lambda =d_p=0$ is also performed in the same way as done for the model (1).

Different timings for taking both the first restriction measure and for lifting it are then investigated, using all the above alternatives.

Finally an intermittent restrictive policy is examined, for which when the infected are observed to trespass a threshold, distancing measures are taken. Here again lock-down or total isolation produce essentially the same results. The use of different timings for the introduction of the restrictions is also scrutinized.