Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan.

We propose a compartmental mathematical model for the spread of the COVID-19 disease with special focus on the transmissibility of super-spreaders individuals. We compute the basic reproduction number threshold, we study the local stability of the disease free equilibrium in terms of the basic reproduction number, and we investigate the sensitivity of the model with respect to the variation of each one of its parameters. Numerical simulations show the suitability of the proposed COVID-19 model for the outbreak that occurred in Wuhan, China.

Introduction

Mathematical models of infectious disease transmission dynamics are now ubiquitous. Such models play an important role in helping to quantify possible infectious disease control and mitigation strategies [1], [2], [3]. There exist a number of models for infectious diseases; as for compartmental models, starting from the very classical SIR model to more complex proposals [4].

Coronavirus disease 2019 (COVID-19) is an infectious disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The disease was first identified December 2019 in Wuhan, the capital of Hubei, China, and has since spread globally, resulting in the ongoing 2020 pandemic outbreak [5]. The COVID-19 pandemic is considered as the biggest global threat worldwide because of thousands of confirmed infections, accompanied by thousands deaths over the world. Notice, by March 26, 2020, report 503,274 confirmed cumulative cases with 22,342 deaths. At the time of this revision, the numbers have increased to 1,353,361 confirmed cumulative cases with 79,235 deaths, according to the report dated by April 8, 2020, by the Word Health Organization.

The global problem of the outbreak has attracted the interest of researchers of different areas, giving rise to a number of proposals to analyze and predict the evolution of the pandemic [6], [7]. Our main contribution is related with considering the class of super-spreaders, which is now appearing in medical journals (see, e.g., [8], [9]). This new class, as added to any compartmental model, implies a number of analysis about disease free equilibrium points, which is also considered in this work.

The manuscript is organized as follows. In Section 2, we propose a new model for COVID-19. A qualitative analysis of the model is investigated in Section 3: in Section 3.1, we compute the basic reproduction number R 0 of the COVID-19 system model; in Section 3.2, we study the local stability of the disease free equilibrium in terms of R 0. The sensitivity of the basic reproduction number R 0 with respect to the parameters of the system model is given in Section 4. The usefulness of our model is then illustrated in Section 5 of numerical simulations, where we use real data from Wuhan. We end with Section 6 of conclusions, discussion, and future research.

The proposed COVID-19 compartment model

Based on a 2016 model [10], and taking into account the existence of super-spreaders in the family of corona virus [11], we propose a new epidemiological compartment model that takes into account the super-spreading phenomenon of some individuals. Moreover, we consider a fatality compartment, related to death due to the virus infection. In doing so, the constant total population size N is subdivided into eight epidemiological classes: susceptible class (S), exposed class (E), symptomatic and infectious class (I), super-spreaders class (P), infectious but asymptomatic class (A), hospitalized (H), recovery class (R), and fatality class (F). The model takes the following form: with β quantifying the human-to-human transmission coefficient per unit time (days) per person, quantifies a high transmission coefficient due to super-spreaders, and l quantifies the relative transmissibility of hospitalized patients. Here κ is the rate at which an individual leaves the exposed class by becoming infectious (symptomatic, super-spreaders or asymptomatic); ρ 1 is the proportion of progression from exposed class E to symptomatic infectious class I; ρ 2 is a relative very low rate at which exposed individuals become super-spreaders while is the progression from exposed to asymptomatic class; γ is the average rate at which symptomatic and super-spreaders individuals become hospitalized; γ is the recovery rate without being hospitalized; γ is the recovery rate of hospitalized patients; and δ, and δ are the disease induced death rates due to infected, super-spreaders, and hospitalized individuals, respectively. At each instant of time, gives the number of death due to the disease. The transmissibility from asymptomatic individuals has been modeled in this way since it was not apparent their behavior. Indeed, at present, this question is a controversial issue for epidemiologists. A flowchart of model (1) is presented in Fig. 1 .

Fig. 1

Flowchart of model (1).

Qualitative analysis of the model

One of the most significant thresholds when studying infectious disease models, which quantifies disease invasion or extinction in a population, is the basic reproduction number [12]. In this section we obtain the basic reproduction number for our model (1) and study the locally asymptotically stability of its disease free equilibrium (see Theorem 1).

The basic reproduction number

The basic reproduction number, as a measure for disease spread in a population, plays an important role in the course and control of an ongoing outbreak. It can be understood as the average number of cases one infected individual generates, over the course of its infectious period, in an otherwise uninfected population. Using the next generation matrix approach outlined in van den Driessche and Watmough [13] to our model (1), the basic reproduction number can be computed by considering the below generation matrices F and V, that is, the Jacobian matrices associated to the rate of appearance of new infections and the net rate out of the corresponding compartments, respectively, where

The basic reproduction number R 0 is obtained as the spectral radius of precisely,

For the parameters used in our simulations (see Table 1), one computes this basic reproduction number to obtain . This means that the epidemic outbreak that has occurred in Wuhan was well controlled by the Chinese authorities.

Table 1

Values of the model parameters corresponding to the situation of Wuhan, as discussed in Section 5, for which .

Name	Description	Value	Units
β	Transmission coefficient from infected individuals	2.55	day^−1
l	Relative transmissibility of hospitalized patients	1.56	dimensionless
β^^′	Transmission coefficient due to super-spreaders	7.65	day^−1
κ	Rate at which exposed become infectious	0.25	day^−1
ρ1	Rate at which exposed people become infected I	0.580	dimensionless
ρ2	Rate at which exposed people become super-spreaders	0.001	dimensionless
γa	Rate of being hospitalized	0.94	day^−1
γi	Recovery rate without being hospitalized	0.27	day^−1
γr	Recovery rate of hospitalized patients	0.5	day^−1
δi	Disease induced death rate due to infected class	3.5	day^−1
δp	Disease induced death rate due to super-spreaders	1	day^−1
δh	Disease induced death rate due to hospitalized class	0.3	day^−1

Local stability in terms of the basic reproduction number

Noting that the two last equations and the fifth of system (1) are uncoupled to the remaining equations of the system, we can easily obtain, by direct integration, the following analytical results:

Furthermore, since the total population size N is constant, one has

Therefore, the local stability of model (1) can be studied through the remaining coupled system of state variables, namely, the variables E, I, P, and H in (1). The Jacobian matrix associated to these variables of (1) is the following one: where ϖ, ϖ, and ϖ are defined in (3). The eigenvalues of the matrix J are the roots of the following characteristic polynomial: where

Next, by using the Liénard–Chipard test [14], [15], all the roots of Z(λ) are negative or have negative real part if, and only if, the following conditions are satisfied:

a > 0, ;

a 1 a 2 > a 3.

In order to check these conditions of the Liénard–Chipard test, we rewrite the coefficients a 1, a 2, a 3, and a 4 of the characteristic polynomial in terms of the basic reproduction number given by (4):

Moreover, we also compute, in terms of R 0, the following expression:

From these previous expressions, it is clear that if R 0 < 1, then the conditions of the Liénard–Chipard test are satisfied and, as a consequence, the disease free equilibrium is stable. In the case when R 0 > 1, we have that a 4 < 0 and, by using Descartes’ rule of signs, we conclude that at least one of the eigenvalues is positive. Therefore, the system is unstable. In conclusion, we have just proved the following result:

The disease free equilibrium of system (1) , that is, (N, 0, 0, 0, 0, 0, 0, 0), is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Next we investigate the sensitiveness of the COVID-19 model (1), with respect to the variation of each one of its parameters, for the endemic threshold (4).

Sensitivity analysis

As we saw in Section 3, the basic reproduction number for the COVID-19 model (1), which we propose in Section 2, is given by (4). The sensitivity analysis for the endemic threshold (4) tells us how important each parameter is to disease transmission. This information is crucial not only for experimental design, but also to data assimilation and reduction of complex models [16]. Sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values, since there are usually errors in collected data and presumed parameter values. It is used to discover parameters that have a high impact on the threshold R 0 and should be targeted by intervention strategies. More accurately, sensitivity indices’ allows us to measure the relative change in a variable when a parameter changes. For that purpose, we use the normalized forward sensitivity index of a variable with respect to a given parameter, which is defined as the ratio of the relative change in the variable to the relative change in the parameter. If such variable is differentiable with respect to the parameter, then the sensitivity index is defined as follows.

See [17], [18]

The normalized forward sensitivity index of R 0, which is differentiable with respect to a given parameter θ, is defined by

The values of the sensitivity indices for the parameters values of Table 1, are presented in Table 2 .

Table 2

Sensitivity of R0 evaluated for the parameter values given in Table 1.

Parameter	Sensitivity index
β	0.963
l	0.631
β^^′	0.366
κ	0.000
ρ1	0.941
ρ2	0.059
γa	0.418
γi	−0.061
γr	−0.395
δi	−0.699
δp	−0.027
δh	−0.238

These values have been determined experimentally in such a way the mathematical model describes well the real data, giving rise to Figs. 2 and 3 . Other values for the parameters can be found, e.g., in Aguilar et al. [19].

Fig. 2

Number of confirmed cases per day. The green line corresponds to the real data obtained from reports [5], [20], [21] while the black line () has been obtained by solving numerically the system of ordinary differential Eq. (1), by using the Matlab code ode45. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Number of confirmed deaths per day. The red line corresponds to the real data obtained from reports [5], [20], [21] while the black line has been obtained by solving numerically, using the Matlab code ode45, our system of ordinary differential Eq. (1) to derive D(t) given in (2). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Note that the sensitivity index may depend on several parameters of the system, but also can be constant, independent of any parameter. For example, means that increasing (decreasing) θ by a given percentage increases (decreases) always R 0 by that same percentage. The estimation of a sensitive parameter should be carefully done, since a small perturbation in such parameter leads to relevant quantitative changes. On the other hand, the estimation of a parameter with a rather small value for the sensitivity index does not require as much attention to estimate, because a small perturbation in that parameter leads to small changes.

From Table 2, we conclude that the most sensitive parameters to the basic reproduction number R 0 of the COVID-19 model (1) are β, ρ 1 and δ. In concrete, an increase of the value of β will increase the basic reproduction number by 96.3% and this happens, in a similar way, for the parameter ρ 1. In contrast, an increase of the value of δ will decrease R 0 by 69.9%.

Numerical simulations: the case study of Wuhan

We perform numerical simulations to compare the results of our model with the real data obtained from several reports published by WHO [20], [21] and worldometer [5].

The starting point of our simulations is 4 January 2020 (day 0), when the Chinese authorities informed about the new virus [20], with already 6 confirmed cases in one day. From this period up to January 19, there is less information about the number of people contracting the disease. Only on January 20, we have the report [21], with 1460 new reported cases in that day and 26 the dead. Thus, the infection gained much more attention from 21 January 2020, with 1739 confirmed cases and 38 the dead, up to 4 March 2020, when the numbers in that day were as low as 11 and 7, respectively infected and dead, after a pick of 3892 confirmed cases on 27 January 2020 and a pick of 254 dead on 4 February 2020. Here we follow the data of the daily reports published by [5]. We show that our COVID-19 model describes well the real data of daily confirmed cases during the 2 months outbreak (66 days to be precise, from January 4 to March 9, 2020).

The total population of Wuhan is about 11 million. During the COVID-19 outbreak, there was a restriction of movements of individuals due to quarantine in the city. As a consequence, there was a limitation on the spread of the disease. In agreement, in our model we consider, as the total population under study, . This denominator has been determined in the first days of the outbreak and later has been proved to be a correct value: according to the real data published by the WHO, it is an appropriate value for the restriction of movements of individuals. As for the initial conditions, the following values have been fixed: and .

We would like to mention that there exist gaps in the reports of the WHO at the beginning of the outbreak. For completeness, we give here the list L of the number of confirmed cases in Wuhan per day, corresponding to the green line of Fig. 2, and the list L of the number of dead individuals in Wuhan per day, corresponding to the red line of Fig. 3:

Lists L and L have 66 numbers, where L(0) represents the number of confirmed cases 04 January 2020 (day 0) and L(65) the number of confirmed cases 09 March 2020 (day 65) and, analogously, L(0) represents the number of dead on January 4, and L(65) the number of dead on March 9, 2020.

Conclusions and discussion

Classical models consider SIR populations. Here we have taken into consideration the super-spreaders (P), hospitalized (H), and fatality class (F), so that its derivative (see formula (2)) gives the number of deaths (D). Our model is an ad hoc compartmental model of the COVID-19, taking into account its particularities, some of them still not well-known, giving a good approximation of the reality of the Wuhan outbreak (see Fig. 2) and predicting a diminishing on the daily number of confirmed cases of the disease. This is in agreement with our computations of the basic reproduction number in Section 4 that, surprisingly, is obtained less than 1. Moreover, it is worth to mention that our model fits also enough well the real data of daily confirmed deaths, as shown in Fig. 3.

Our theoretical findings and numerical results adapt well to the real data and it reflects or reflected the reality in Wuhan, China. The number of hospitalized persons is relevant to give an estimate of the Intensive Care Units (ICU) needed. Some preliminary simulations indicate that this would be useful for the health authorities. Our model can also be used to study the reality of other countries, whose outbreaks are currently on the rise. We claim that some mathematical models like the one we have proposed here will contribute to reveal some important aspects of this pandemia.

Of course, this investigation has some limitations, being the first on the relative recent spread of the new coronavirus and therefore the limited data accessible at the beginning of this study. In the future, we can develop further this prototype. Even with these shortcomings, the model can be useful due to the high relevance of the topic. Finally, we suggest new directions for further research:

the transmissibility from asymptomatic individuals;

compare, in the near future, our results with other models;

consider sub-populations related to age, gender, etc.;

introduce preventive measures in this COVID-19 epidemic and for future viruses;

integrate into the model some imprecise data by using fuzzy differential equations;

include the viral load of the infectious into the model.

These and other questions are under current investigation and will be addressed elsewhere.

Funding

This research was funded by the (FCT) within project UIDB/04106/2020 (CIDMA). Ndaïrou is also grateful to the support of FCT through the Ph.D. fellowship PD/BD/150273/2019. The work of Area and Nieto has been partially supported by the (AEI) of Spain, cofinanced by the (FEDER) corresponding to the 2014–2020 multiyear financial framework, project MTM2016-75140-P. Moreover, Nieto also thanks partial financial support by under grant ED431C 2019/02.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.