A simple but complex enough θ -SIR type model to be used with COVID-19 real data. Application to the case of Italy.

Since the start of the COVID-19 pandemic in China many models have appeared in the literature, trying to simulate its dynamics. Focusing on modeling the biological and sociological mechanisms which influence the disease spread, the basic reference example is the SIR model. However, it is too simple to be able to model those mechanisms (including the three main types of control measures: social distancing, contact tracing and health system measures) to fit real data and to simulate possible future scenarios. A question, then, arises: how much and how do we need to complexify a SIR model? We develop a θ -SEIHQRD model, which may be the simplest one satisfying the mentioned requirements for arbitrary territories and can be simplified in particular cases. We show its very good performance in the Italian case and study different future scenarios.

Introduction

The outbreak of novel coronavirus disease 2019 (COVID-19) continues to spread rapidly around the world [1]. Originated in Wuhan in central China, where it was already active at the beginning of December 2019, reached multiple countries in merely a month until it became a pandemic as of 11 March 2020 [2] when the disease was confirmed in more than 118,000 cases reported globally in 114 countries.

Although some authors (see, e.g., [3]) propose SEIR type models with little variations, COVID-19 is a disease caused by a new virus and needs a model taking into account its known specific characteristics. In particular, it would be convenient to develop a model which satisfies the following: (a) includes the biological and sociological mechanisms influencing the disease spread, (b) its parameters have a clear interpretation in terms of those mechanisms (such as a contact rate, fatality ratios, duration of infection, etc.), and (c) incorporates

the effect of undetected infected people (see [4]), being able to show the dependence of the impact of COVID-19 on the percentage of detected cases over the real total infected cases,

the effect of different infectiousness conditions of hospitalized people (differentiating those with mild and severe conditions that will recover from those who will finally die),

the estimation of the needs of beds in hospitals (which is one of the major problems for policy makers addressing COVID-19),

the effect of different control measures.

Some of the models that appeared at the beginning of the pandemic, especially for the study of the first wave (see [3], [4], [5], [6]), do not incorporate all the special characteristics (undetected cases, undetected deaths, quarantine, social distancing, etc.) of the COVID-19 pandemic. Some of them are very complex in terms of the number of compartments and connections between them, trying to simulate many details that could probably be neglected. Furthermore, after developing a model, it would be desirable to check how it fits real data and to use it in order to study future scenarios with different control measures (here, we consider social distancing, contact tracing and health system measures).

This paper does not take into account the possibility of resusceptibility as it is done, for example, in [7]. This is due to the fact that the number of confirmed reinfections is still not relevant enough compared with the total figures of the pandemic. According to [8], only 27 cases of reinfection were documented worldwide when it was written. Furthermore, the webpage https://bnonews.com/index.php/2020/08/covid-19-reinfection-tracker/ reports an updated list of reinfections, with a total of 31 cases, at the time we wrote this article.

In [9], a disease contact rate () and with different values before and after lockdown of countries are considered and applied in a SEIR type model for the confirmed cases tested with several countries. The model presented here allows to model not only confirmed cases but also deaths, hospitalized and recovered by considering different values of that depend on the particular control measures for every time period.

Concerning the particular case of Italy, some studies have dealt with its particular situation, such as [10], [11], [12] and [13]. It is important to note that all the models presented in these works only serve the particular case of Italy, and are not general models applicable to other countries.

The paper [10] is an initial work published in mid-March 2020, which highlights the importance of this pandemic and the serious consequences it can have in countries like Italy where data show rapid progress in infections, which can cause a health system collapse due to the lack of beds in hospitals and ICU posts. The data are fitted using an exponential type model, which is not realistic for long-term modeling, although it allows a wake-up call to the country’s health authorities.

In [11], the distribution in terms of mobility is also studied through a network of the different Italian regions. It is confirmed the effectiveness at a regional level of the national lockdown strategy and proposed coordinated regional interventions to prevent future national lockdowns. However, asymptomatic patients are not considered.

In [12], a geographical model is presented to describe the number of daily cases, hospitalized and deaths within the Italian regions for the first month of the disease and study the effects of drastic measures for transmission containment. Authors concentrate on estimating the effects of severe progressive restrictions posed to human mobility and human-to-human contacts in Italy. As a novelty, a presymptom compartment is considered. The proposed model is extended and further developed in [13], where the effect of mobility between regions and provinces within Italy is also considered, taking into account the control measures at all times to describe the variable of new hospitalizations. The number of people who have to be isolated in each region is calculated to avoid a drastic growth in the number of infected. This model is very interesting, because it incorporates the variable of mobility between regions, but only the number of new hospitalizations is described, without studying the number of deaths since they focus on the implications of the relaxation of confinement measures.

An interesting model appearing in [5], tested for the case of Italy, has some similarities with the one presented here, but there are some important differences to be remarked. First of all, the model we propose is a general one that may be used with any country or region, just considering the correspondent data, whereas the model in [5] seems to be defined for the particular case of the first wave in Italy, although it could be adapted to other territories. Secondly, our model considers the compartment for people in quarantine (, detected asymptomatic or with mild symptoms, infectious people) that is not the same as the proposed compartment for diagnosed (, asymptomatic infected, detected) in the former model, which is even connected to the possibility of dying. Moreover, the general role of the undetected people in our model is completely different and some of them can even die (those who are in ).

In this paper, we develop a general deterministic epidemiological model able to simulate the biological and sociological COVID-19 mechanisms in different territories, to fit real data, and to simulate different future scenarios of control measures. The choice of using a deterministic model instead of a stochastic one is done as a first approach, since such kind of models present some advantages, such as: a low computational complexity allowing a better calibration of the model parameters or the possibility of using the theory of ordinary differential equations for suitably analyzing and interpreting the model. Furthermore, according to [14], deterministic models should be the first tool to be used when modeling a new problem. The authors of that work also note that the stochastic models are not suitable when it is difficult or impossible to determine the probability distribution, are difficult to analyze and require more data for the calibration of the model. We note that, although we use a deterministic model, we must deal with some probabilistic processes. In those cases, random values are substituted by their expected values.

We look for the simplest model, having the minimum complexity to take into account all the main non-negligible processes of this pandemic. When simulating the dynamics of the COVID-19 disease in a territory, our model can be used in its full general form or in simplified particular versions of it, depending on the way the territory deals with the pandemic. Of course, the model may be modified by introducing age structure, stochasticity, levels of uncertainty, mobility between territories, spatial structure (see, e.g., [15]), etc. in order to focus in special features of the pandemic or to refine some results.

The reference epidemiological model is the SIR one, which divides the population of a territory in three compartments: Susceptible (), Infectious () and Recovered (). For the COVID-19 pandemic these three compartments are not enough. First of all, the incubation period for the new coronavirus SARS-CoV-2 is between 5 and 6 days (see [16]), which is important to incorporate in the model by means of a compartment for Exposed () people. Secondly, estimation of needs of beds in hospital is important for policy makers, which implies the convenience of suitable compartments for people in Hospitals (), following a compartment () for an initial stage of the infectious period. Moreover, the disease contact rates inside and outside hospitals are usually different. Thirdly, one of the main data reported each day is the number of Deaths (), which deserves another compartment and makes convenient to split in two branches, one for people in hospital that will recover () and the other one for people in hospital that will die (). Therefore, we need, at least, a SEIHRD model.

Furthermore, one of the characteristics of COVID-19 is the existence of undetected infectious cases, either because they are asymptomatic (see [4]) or because, even if they have symptoms, the corresponding country/territory does not have enough tests to check the infectiousness or the health system is saturated. In order to be able to consider the impact of this COVID-19 characteristic, two different branches should distinguish between detected and undetected cases and deaths. This can be done in terms of the instantaneous ratio of new detected infectious cases per unit time over the total (detected and undetected) new infectious cases per unit time, generating the corresponding -SEIHRD model. As done previously with people in hospital, a distinction can be done between the undetected infectious people (compartment ) that will recover (compartment ) and the undetected infectious people (compartment ) that will die (compartment ). Finally, some territories (like Italy) use a quarantine period (typically at home) for people leaving the hospital and for detected asymptomatic (or with mild symptoms) infectious people, which makes necessary a new compartment (Q), which completes the full -SEIHQRD general model. A diagram summarizing the model is given in Fig. 1.

Fig. 1

Diagram summarizing the model for COVID-19 given by system (1).

A detailed description and discussion of the model, its parameters and outputs is provided in Sections 2, 3 . We point out that all the parameters in the model have well-understood epidemiological interpretations such as average duration of incubation period, contact rates, etc. In Section 2, we also show some of the outputs of the model and how to compute them. Particular outputs of interest for epidemiologists are the computation of the associated basic reproduction number and the evolution of the effective reproduction number.

Once we have developed the model, the question is how to use it with real cases. Let us suppose we want to use it to study the dynamics of the COVID-19 in a particular territory. Three steps can be followed:

The first thing is to decide if all the compartments are necessary, or only some of them. For instance, in China and Senegal all infected people are treated in hospitals (or other kind of health facilities) and, therefore, compartment is not necessary [6]. However, other territories, like Italy or Spain, use the quarantine at home for infected people with very mild symptoms or not symptoms at all. Other possibilities of simplifying the model are also possible, for particular territories.

After deciding the compartments that will be used with a particular territory, the following step is to choose the parameters involved in the model. To do that we use values given in the literature (as the duration of the incubation period), parametrization of functions involved in the model and a multiobjective technique, to identify parameters. In Section 3 we explain how to do it. This is, actually, part of the validation process, consisting in comparing the results of the simulations (run with the chosen parameters) with the available data about the impact of the pandemic in the territory under study. Typically, a good fitting of simulations outputs to official temporal series (of cases, deaths, people in hospital, etc.) is necessary.

Finally, if the validation process is successful, we can simulate different future scenarios to study the efficiency of some control measures (being increased or relaxed), the impact of the undetected cases in the dynamics of the pandemic, etc.

In Section 4, we illustrate those three steps with the particular case of Italy.

Mathematical formulation of the model

In this section we give a detailed description of the -SEIHQRD model that we propose for the COVID-19. Then, we described some outputs of the model that will be used in the numerical experiments performed in Section 4.

The mathematical model

According to the known characteristics of the COVID-19 pandemic, we assume that each person is in one of the following compartments (see [17], [18], [19], [20]):

Susceptible (denoted by ): The person is not infected by the disease pathogen.

Exposed (denoted by ): The person is in the incubation period after being infected by the disease pathogen, and has no clinical signs. This stage is also called presymptomatic period. According to the available literature (see, e.g., [17]), there exist some evidences that a presymptomatic person could infect other people but with a lower probability than people in the infectious compartments. After the incubation period, the person passes to the Infectious compartment .

Infectious (denoted by ): After the incubation period, it is the first compartment of the infectious period, where nobody is expected to be detected yet. The person has finished the incubation period and may infect other people (even if he/she is asymptomatic, see, e.g., [19]). After this period, people in this compartment can be, either taken in charge by health authorities (and we classify them as hospitalized or quarantine, according to the severity of their symptoms), or not detected by authorities and continue as infectious (but in another compartment, or ).

Infectious but undetected that will survive the disease (denoted by ): After being in the compartment , the person can still infect other people but is not detected and is not reported by authorities. We assume that only asymptomatic infected people or those with low or medium symptoms can reach this compartment, not the people who will die. After this period, people in this compartment survive the disease and pass to the Recovered compartment .

Infectious but undetected that will die (denoted by ): After being in the compartment , the person can still infect other people, have clinical signs but not be detected and reported by authorities. We assume that a person in this state is part of a population that is highly vulnerable to the COVID-19 (such as people in retirement homes) and will die quickly without being treated in hospital. After being in this state, a person passes to the compartment .

We remark that vulnerable people died outside hospitals during the first wave of the pandemic in some countries (such as Italy or Spain; see [21], [22], [23], [24], [25]). Indeed, due to the saturation of the whole health system (mainly because of lack of beds and/or healthcare workers), this vulnerable population did not receive adequate medical cares and quickly died before being hospitalized. In some of those countries (such as Italy, see [21], [22]), a few a-posteriori estimations of the number of people dead due to COVID-19 and not detected have been reported by authorities or other organizations and, currently, efforts are made to avoid reproducing this extreme situation (see [25]). Regarding those few data and for the sake of simplicity of the model, we have decided to model people who have suffered this situation by using the specific compartment . We note that the results reported later in Section 4 indicate that this choice is adequate.

Quarantine (denoted by ): The person is in quarantine at home or some special facilities (different from hospitals) during a fixed period.

For example, according to [26], [27], the isolation is performed in one room, always wearing mask at home, etc. to reduce the risk of transmission to cohabiting people. Of course, infection could occur between people in the same home, but this is usually before entering the quarantine period (i.e., while being in the compartment ) and the model takes that into account. Furthermore, according to [5], transmission between people in quarantine at home (thus, aware of their infection) and other person (including cohabiting persons) is severely limited.

Therefore, in this work we assume that a person in compartment cannot infect other people (although other possibilities are also easy to include, as explained below) and at the end of this state passes to the compartment . The use of this compartment in the model depends on the strategy of each territory (some of them do not use this control strategy).

Hospitalized that will recover (denoted by ): The person is in hospital and (according to [5], [28], [29]) can still infect other people (in particular healthcare workers). At the end of this state, a person passes to the compartment (if implemented in the territory under study; otherwise the person passes directly to the compartment ).

Hospitalized that will die (denoted by ): The person is hospitalized and can still infect other people. At the end of this state, a person passes to the compartment .

Dead by COVID-19 (denoted by ) after being at hospital.

Dead by COVID-19 but undetected (denoted by ): The person has not survived the disease and the death was not classified by authorities as caused by COVID-19.

Recovered after being previously detected as infectious (denoted by ): The person was previously detected as infectious, survived the disease, passed the quarantine period (if implemented in the territory under study), is no longer infectious and has developed a natural immunity to the virus. Although some cases of reinfection are known, following [18], they seem to be isolated cases; thus, despite the fact that the evolution of this immunity is unknown for a longer period of time, we make the natural immunity assumption. However, it can be modified if new findings regarding this topic are found.

Recovered after being previously infectious but undetected (denoted by ): The person was not previously detected as infectious, survived the disease, is no longer infectious and has developed a natural immunity to the virus (with similar assumptions as above).

We remark that, in this work, we work with two compartments, namely and described previously, to model the dynamics of the number of people in hospital. One may note that a more detailed classification of hospitalization could be considered (for instance, according to the severity of the symptoms, the availability of medical material or beds, etc.) to give more recommendations to decision makers. This could be done by introducing additional compartments to the model. However, this would require the use of detailed data regarding people in hospital which, to the best of our knowledge, are not publicly available for most of the countries (in particular, Italy). Furthermore, adding more compartments will complexify deeply the model, making difficult the parameter estimation process. For those reasons, we decided to use two compartments. We highlight the fact that the results reported in Section 4 (see Fig. 2) indicate that this choice is reasonable as the model estimates quite well the reported data related to people in hospital (which are, for the case of Italy, the number of people in hospital, without any specific classification).

Fig. 2

Comparison of some of the outputs of the simulation run with (Est.) and the Italian official reported data (Rep.). Top-left: Cumulative number of cases and deaths. Top-right: People in hospital. Bottom-left: People in quarantine. Bottom-right: Cumulative number of recovered people.

The authorities may apply various control measures in order to control the COVID-19 spread (see [30]):

Isolation: Infected people are isolated from contact with other people. Only health professionals are in contact with them. However, infection of those professionals also occurs (see [16]). Isolated patients receive an adequate medical treatment that reduces the COVID-19 fatality rate.

Lockdown areas: Movement of people in areas with a high number of infected people may be restricted and controlled, to avoid the spread of the disease.

Contact tracing: The objective of tracing is to identify potential infectious contacts which may have infected a person or spread COVID-19 to other people. Increase the number of tests in order to increase the percentage of detected infected people.

Increase of health system resources: The number of operational beds and personal available to detect and treat affected people is increased, producing a decrease in the length of the infectious period for the compartment .

The model is used to evaluate the spread of a human disease within a considered territory during a fixed time interval.

At the beginning of the simulation, the model parameters are set by the user (as shown below). We can start our simulation at any initial time , considering that only susceptible people live in the territories that are free of the disease, whereas the number of people in compartments , , , , , , , , , and of the infected territories are set to their corresponding values. Then, for the time interval , with being the maximum number of simulation days, the model is applied. If at the end of a simulation day there is less than 1 person in each compartment , , , , and in all the territory, the simulation is stopped. Else, the simulation is stopped when . Furthermore, the control measures are also implemented and they can be activated or deactivated, when starting the model, in order to quantify their effectiveness to reduce the magnitude and duration of the COVID-19 pandemic in the territory under study.

A diagram summarizing the main structure of our full model is presented in Fig. 1.

We assume that people in a territory are characterized to be in one of the compartments described above: , , , , , , , , , , or . For the sake of simplicity, we assume that, at each time, the population inside a territory is homogeneously distributed. Thus, the spatial distribution of the epidemic inside a territory is omitted (it can be taken into account by dividing some territories into a set of smaller regions with similar characteristics). Since the natural natality and mortality (not from COVID-19) do not seem to be important factors for COVID-19 (at least for relatively short periods of time), they are omitted in the model (actually, they can be readily included, if necessary).

In this work we will only consider the within-country/territory disease spread of one territory where, starting from suitable values of , the COVID-19 pandemic is already spreading by its own, with a relatively negligible dependence on the movement of people in/out of this territory. The between-country/territory disease spread may be also taken into account in the model following the approach developed in [15] with the Be-CoDiS model.

There are some studies about the relevance of humidity and temperature for the spread of COVID-19. In [31], it is shown that the observed patterns of COVID-19 are not completely consistent with the hypothesis that high absolute humidity may limit the survival and transmission of this virus. Furthermore, in [32] it has been found that the lower is the temperature, the greater is the survival period of the SARS-CoV-2 outside the host. Since there is no clear scientific evidence of the effect of the humidity and the temperature on the SARS-CoV-2, we have not included these two factors in our model (this would need to be revised in case of appearing new findings regarding this topic).

Under those assumptions, the evolution of the compartments mentioned above is modeled by the following system of ordinary differential equations: In this system we have used, as stated above, that people in compartment cannot infect other people. To include the possibility of infection, we would only need to add the term inside the parenthesis of the first two equations.

In (1):

is the number of people in the territory before the start of the pandemic.

is the instantaneous infection detected fatality ratio (iIdFR) in the territory, at time : the proportion of new detected people (after being in compartment ) that will die of COVID-19 after days, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time. More precisely, where (respectively, ) is the model cumulative number of COVID-19 detected (respectively, undetected) cases at day (see more details about them below, when describing the outputs of the model).

is the instantaneous infection undetected fatality ratio (iIuFR) in the territory, at time : the proportion of new undetected people (after being in compartment ) that will die of COVID-19 after days, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time. More precisely,

We point out that is the instantaneous infection fatality ratio (iIFR), at time : the proportion of new (detected or undetected) people (after being in compartment I) that will die of COVID-19, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time. This value is closely related to the IFR defined below (see the Results section). Actually, if it is assumed that the number of total deaths is proportional to the number of total cases, then IFR is, essentially, .

is the proportion of the number of new infectious people that are detected and documented by the authorities (after being in compartment ), per unit time, compared to the total number of new infectious people (detected or undetected), per unit time, at time . More precisely, We point out that is the instantaneous case fatality ratio (iCFR), at time : the proportion of the number of new detected people (after being in compartment ) that will die of COVID-19 after days, per unit time, compared to the number of new detected people (after being in compartment I), per unit time.

This value is closely related to the Case Fatality Ratio (CFR) defined as the expected value of the probability of dying for a person who is infected and detected, i.e. the number of detected deaths due to the virus over the number of detected infections. Actually, if it is assumed that the number of detected deaths is proportional to the number of detected cases, then CFR is, essentially, .

are the disease contact rates (day−1) of a person in the respective compartments. They may change with time, due to the control measures.

is the transition rate (day−1) from compartment to compartments , , , or , at time . Similarly, , and denote the transition rates (day−1) from compartments , , , , and to compartments , , , , and , respectively, at time . For some particular territories, some compartments could be excluded from the model.

is the ratio of the number of new detected infected people that will survive the disease and are hospitalized (they go to compartment ) over the total number of new detected infected people that will survive the disease (and may also do quarantine, compartment , without going to hospital).

and are the people infected that arrive/leave from/to other territories, per day. Both can be modeled following the between-country spread part of the Be-CoDiS model described in [15]. For the sake of simplicity, we have included those terms only in the second equation, corresponding to compartment , but they could be included in the equations for other compartments too. When a discrete sequence of imported/exported infections is known it can be included in the numerical resolution of the system by simple addition/removal to/from the suitable compartment, at the corresponding discrete time.

System (1) is completed with initial data , , , , , , , , , , and given in . These values must be non-negative and . It is standard (see, e.g., [33], [34]) to prove that , , , , , , , , , , and are non-negative and , for all .

We remark that system (1) is an epidemic-type model (as described, e.g., in [35]), which is not trying to capture endemic behaviors. Actually, there is no incoming flux in the Susceptible compartment; therefore, the fate of every individual will be either staying susceptible, recovering or dying. This means there are no endemic equilibria in our system; on the contrary, there exists a manifold containing infinite disease-free equilibrium points.

We point out that the 8th, 9th, 10th, 11th and 12th equations of system (1) are not coupled with the other equations. Thus, we can solve the first seven equations of that system and the solution of the last five ones can be computed as follows: is the solution of the 8th equation of system (1), which is a linear ordinary differential equation (we do not include here the easy corresponding explicit solution, for the sake of simplicity) and It is important to remark that, if for some territory the model does not include the compartment , then , and

For the numerical simulations presented in this paper, the first eight equations of system (1) were numerically solved with the classic 4th order Runge–Kutta method with four stages (RK4) and with four hours as time step (which was tested to see that it is suitable to get stable results).

Outputs of the model

Here, we present some of the outputs that can be obtained from the mathematical model, used to analyze the results of the simulations:

: The model cumulative number of COVID-19 detected cases at day , which is given by and can be also computed as follows:

: The number of new cases per day at day , which is given by

: The model cumulative number of COVID-19 detected infected healthcare workers at day , which is given by

: The model cumulative number of COVID-19 undetected cases at day , which is given by and can be also computed as follows:

: The model cumulative number of detected deaths (due to COVID-19), at day , which is given by .

: The model cumulative number of undetected deaths, at day , which is given by .

R: The effective reproduction number of COVID-19. It is defined as the number of cases that one infected person generates, on average, over the course of its infectious period [36]. Part of the population can be already infected and/or control measures may have been implemented. It depends on the considered territory and changes during the spread of the disease. Furthermore, if at time no control measures have been applied, R R, where is the basic reproduction number of COVID-19, which is defined as the number of cases that one infected person generates on average over the course of its infectious period, in an otherwise uninfected population and without special control measures. Typically, the spread of the disease slows down when R. Taking into account this definition, for cases with constant transition rates we propose the following non-classical estimation: In this formula, we estimate the mean number of secondary cases produced by one person infected at time , during its whole infectiousness period, taking into account all the possible evolutions of its health status.

Composite rectangular, trapezoidal or Simpson rules can be used to obtain an estimation of (2).

We remark that other estimations for R can be also considered. For instance, one can use the so-called “Next generation Matrix” method to approximate its value (see [37]). However, as pointed in [38], this approximation does not correspond accurately to the biological definition of R (e.g., parameters are assumed to be constant and the change of the proportion of susceptible people during the infectious time of a person is omitted). We remark that R should be used with caution by policymakers, without relying only on that value (see [39]).

Note that we can estimate the particular contribution of each infectious compartment (i.e. , , , , and ) to the value of R and R by computing only their respective integrals in (2).

Hos: The number of people in hospital is estimated as follows: This function can help to estimate and plan the number of clinical beds needed to treat all the COVID-19 cases.

: The number of new hospitalized people per day at day , which is given by

MHos: The maximum number of hospitalized people at the same time (in territory ) during the time interval . It is computed as: This number can help to estimate and plan the number of clinical beds needed to treat all the COVID-19 cases.

and : The number of people infected during the time interval , by contact with people in compartment and , respectively. They are given by:

We point out that those outputs and others that can be computed similarly can be compared with reported data when available.

Model parameter estimation for COVID-19

Here, we present the methodology we used to estimate the parameters of the model described in Section 2.1. We detail each kind of parameter according to its category.

Parametrization of some of the functions involved in the model

The model involves several functions. We now detail how some of these functions can be written, using empirical assumptions, in terms of constant parameters, which will have to be later identified for each territory under study (we will do it below for the case of Italy). Of course, other alternatives may be also used, taking into account the particularities of each territory.

– Disease contact rates: We consider that Here are the respective disease contact rates when no control measures are taken and are functions simulating de control measures, as described below.

Regarding , we may consider where is a fatality rate value that is explained below and and are suitable lower and upper bounds, respectively. For the sake of simplicity, here we take .

Furthermore, following the idea proposed in [4], [6], [15], we assume that there exists a relationship between the contact rates , , and .

More precisely, following [6], we consider that people in compartments , and are less infectious than people in compartment (due to their probably lower virus load or isolation measures, see [4], [17], [19]). Thus, we consider that and and , where .

A similar approach can be followed for and assuming, for instance, that . Alternatively (as we do here for the case of Italy) we can compute good estimates for and using historical data, as shown below.

–Control measures: We may distinguish three different types of control measures: Social distancing, contact tracing and health system measures. Each of them may be included in the model by using different functions, as explained below.

Functions , , (%) represent the reduction of some of the disease contact rates, at time , due to the efficiency of the control measures corresponding to social distancing applied to the compartments , and , respectively (see, e.g., [40]). Here, for the sake of simplicity, we consider that those three functions are equal (of course, other alternatives are also possible) and given by (see [6])

We may also define in a similar way functions and related with contact tracing and health system control measures, respectively, which include similar parameters denoted by and . Here, are the number of changes of control strategies.

Here, for every , and , constants measure the intensity of the control measures (lower value implies higher intensity), (day−1) simulate the efficiency of the control measures (higher value implies faster effectiveness of the control measures) and denote days of changes of each control strategy.

The values of are typically taken from the literature (using the dates when the territories implement special control measures). Some of the values of can be also known. The rest of the parameters need to be identified, as explained below. When doing this, to ease the search process, we use new parameters instead of , as follows:

If we want , then we look for a suitable parameter such that .

If we want , then we look for a suitable parameter such that .

A similar approach can be followed by using new parameters instead of .

Function is related with the contact tracing measures, since the higher those measures are at time , the closer is to 1. The other function involved with the contact tracing measures is , which is used below for estimating the transition rates. The function is used to compute the fatality ratios and to define , representing the reduction of the disease contact rates , at time , due to the efficiency of the health system measures applied to the corresponding compartments. We note that this approach is not taken into account when using the alternative mentioned above to compute estimates for and using historical data (which will be shown below).

–Fatality ratios: As observed in other epidemics (see, e.g. [15]), can be affected by the application of health system measures. Thus, following [6], we propose to consider where is the iIdFR when no health system measures are applied (i.e. ); and is the iIdFR when the implemented health system measures are fully applied (i.e. ). We note that , thus, we consider , with . Alternatively (as we do here for the case of Italy), these constants can be estimated from values of the IFR in the literature (see, e.g., [41], [42]).

Additionally, we assume that where and the corresponding values , are parameters to be chosen for each territory, depending on its particular circumstances.

–Transition rates:

According to [4], the transition rate from to depends only on the disease and, therefore, is considered constant.

Some other transition rates could depend on time. For instance, according to [4], [43], depending on the studied territory, the value of can be increased due to the application of contact tracing control measures (i.e. people with symptoms are detected earlier). As a consequence, the values of , , , and can be decreased (e.g. people with symptoms stay under observation during more time). In order to parametrize these functions, we denote by , , , , and the initial duration (given as a suitable average) in days of a person in compartment , , , , and , respectively, without the application of special contact tracing measures. Then, following [15], we may use the following functions: where represents the decrease of the duration of due to the application of contact tracing measures, at time ; and is the maximum number of days that can be decreased.

We remark that, for some territories, it can be suitable to consider constant (i.e. not dependent on time) all the transition rates, which corresponds to the case with .

–Ratio : We describe it by using the following function depending on the value of (we remark that other alternatives are also possible): Here and , where and are described below (when we show how to estimate ). Indeed, we assume that, if all new detected infected people, when comparing to the case of and , will go to state , as they should present light symptoms. Furthermore, we use a linear interpolation between the case (with , i.e., only persons with severe symptoms are detected) and (with ). More complex formulas could be also considered, but we use the previous one, for the sake of simplicity.

We remark that one of the functions involved in the model, namely , has not been parameterized above. We will show below how to deal with it.

Estimation of parameters and non-parameterized functions involved in the model

Once the model is built, if we want to use it for a particular country or territory, we need to identify suitable model parameters and functions for the case under study. Some of them can be found in the literature. However, due to lack of information regarding the behavior of the SARS-CoV-2 and implicit uncertainty, others need to be estimated.

To do that, we need as much information as possible about the dynamics of the pandemic in the studied territory. In particular, we may use the following temporal series of the territory under study, if available:

: The reported cumulative number of COVID-19 infected, cases at time .

: The reported cumulative number of COVID-19 deaths, at time .

: The reported cumulative number of COVID-19 infected healthcare workers, at time .

: The reported number of known imported cases on day . More precisely, infected people arriving from an external territory. Usually they are in the incubation period, but, in some cases, they are moving between hospitals of different countries or territories (this last possibility is not included in the equations, for the sake of simplicity, but it can be easily done).

: The reported number of known exported cases on day . More precisely, infected people going to an external territory. Usually, they are in the incubation period, but sometimes they are moving between hospitals of different countries or territories (this last possibility is not included in the equations, for the sake of simplicity, but it can be easily done).

: The reported cumulative number of estimated undetected COVID-19 deaths, at time .

: The reported number of people in hospital (i.e. ), at time .

: The number of people in quarantine at home, due to their infection with COVID-19.

We note that, sometimes, there is only partial information and the gaps can be filled with suitable techniques (see, e.g., [44]).

Now, with that information we may estimate some of the model parameters and functions, as follows:

–Ratio of new infected people that are detected: To estimate , first, we compute an estimation of iCFR () as follows: where days is the first date for which and . We note that, in this formula, we estimate the iCFR by taking into account the delay between the detection of a case and its death (instead of the ratio between cases and deaths reported at the same date). This delay could also be considered in the system of differential equations, instead of the exponential-decay-formulation, leading to a system of delay differential equations (see, e.g., [45]).

After that, to regularize the obtained iCFR curve, we apply a filter that consists in computing the average value of the last week. We denote the resulting filtered curve by : To avoid unrealistic situations (due to poor data), if (which is assumed too low, see, e.g., [42]) for some , then we extend the average of the last week used to compute this by including the previous days until getting .

Finally, we define by: where and days to guarantee that is computed by using the average value defined by the second line of (3), for the sake of smoothness. Notice that, if we simulate until the current date, we do not have enough data to compute . In this case, we set to a constant value for the last days of the simulation without data. When simulating scenarios for future dates, we can also keep the last value of or use alternative possible values.

For particular territories (as we do here for the case of Italy) we can use different values of IFR found in the literature and use them as in the previous formula. That will provide different estimations for .

–Contact rates and : We use the following estimation of the fraction of infected healthcare workers (assumed to be due to contact with hospitalized people, see, e.g., [29]): To compute , we first define as and is the first date for which is reported. We point out that, in the previous formula, we have made a distinction about the knowledge of the values of , because they are rarely given every day. We remark that we are taking into account the delay of days between the infection of a person and its detection by authorities. Thus, we set to a constant value for the last days. Finally, to avoid unrealistic situations (due to poor data), we take the following seven-day average value:

Assuming, for the sake of simplicity, that , we have that In the previous expression, the coefficients depend on time .

–Other unknown parameters: We can identify them with a multiobjective technique, as described below.

Identification of some parameters with a multiobjective optimization technique

First, we determine the set of unknown model parameters to be identified by the approach presented here (which will depend on the available information in the considered territory). For each of those parameters, we establish a suitable interval of possible values (given by upper and lower bounds of each parameter) and denote by the feasible region given by the cartesian product of all those intervals, where is the number of unknown parameters.

Once is defined, we follow the methodology proposed in [46]. This approach aims at minimizing the differences between the temporal series reported by the authorities and the corresponding ones obtained as model outputs when considering particular values . Then, given two weight coefficients and and a temporal series , we define the function that measures the error between the reported data, which is denoted by , and the model output denoted by . Notice that, as done in previous works [6], [46], it is based on the relative error. However, as a novelty, we add two new terms to penalize those obtaining a temporal series that underestimate the reported data. We recall that, to reduce the risk of over-fitting to the reported data [47], we use the norm in , which is .

The number of functions depends on the number of available temporal series of the territory under study. The authorities usually report at least the cumulative daily number of cases and deaths (i.e.  and ), see, e.g., [48]. In some cases, they also include other data as, for instance, the number of people in hospital (i.e. ), in quarantine at home (i.e. ), or the cumulative number of undetected deaths (i.e. ) (see, e.g., [49]).

Then, we formulate a multiobjective problem that may have as many objective functions as available temporal series. For instance, we could solve where is the feasible region (i.e. upper and lower bounds of each parameter) and is the number of parameters in . We can also add together (with suitable weights) several functions, in order to reduce the number of functions to be minimized (and the complexity of the corresponding multiobjective problem).

Next, to solve this multiobjective optimization problem, the algorithm called Weighting Achievement Scalarizing Function Genetic Algorithm (WASF-GA) is applied [44], [50], [51]. We choose this particular multiobjective algorithm as it has been applied with success in a previous work to identify some parameters of an epidemiological model simulating the evolution of Ebola Virus Disease outbreaks in Democratic Republic of the Congo [44].

From a general point of view, WASF-GA is an evolutionary algorithm that manages a population of possible solutions (called individuals), trying to improve them iteratively towards the optimal set (called Pareto optimal set). To do that, it uses two procedures at each iteration: reproduction and replacement. Firstly, the reproduction method applies some classical genetic operators (here: selection, mutation, and crossover) to the population of the previous iteration to obtain a new offspring population. Secondly, the replacement procedure decides which solutions from the previous iteration and the offspring populations should compose the population for the next iteration. The main feature of WASF-GA is that this decision is based on preferences, which are defined as reference values for the objectives. Indeed, it classifies the individuals using an achievement scalarizing function (see [51]) characterized by the L distance between the objective function values and the reference values with different weights. The considered sample of weight vectors is built following [51] to ensure they are well distributed.

As a result, WASF-GA returns a set of solutions covering a region of the Pareto optimal front.

In this work, we set the number of iterations to 15 000, the population size to 50 individuals, the number of weight vectors to 50, and the reference values for the objectives to zero, corresponding to the ideal situation of a total matching between the model outputs and the reported data. The remaining input parameters of WASF-GA are those reported in [51]. Those values were set experimentally in order to obtain a good ratio between computational time and precision of the results.

Since WASF-GA is a stochastic algorithm, it may return different Pareto front approximations for different executions of the same instance. Then, to guarantee that results do not depend on the stochasticity of the algorithm, we perform 60 runs of each experiment. Among all the 60 sets of solutions, we extract the most balanced solution, meaning the one whose objective values are the closest to the reference values (in our case, zero) when considering the Euclidean distance.

In the 8th, 9th, 10th, 11th and 12th equations of system (1), which, as said previously, are not coupled with the other equations, there is a parameter, namely , that only appears in those equations. Therefore, its estimation can be done once the multiobjective process is over. For instance, if the authorities report the time series corresponding to the number of people in quarantine at home (i.e. ), the parameter can be estimated separately as a final step of the identification methodology. To do that, we assume that can be parameterized by parameters denoted by (see next section for a particular example). Then, we formulate the following mono-objective problem: where measures the error (as defined by (5)) between the reported data and , the number of people in quarantine returned by the model when considering described by .

Here, we solve this mono-objective problem by using a genetic algorithm configured with the same selection, mutation, and crossover operators as WASF-GA and the same input parameters (i.e. we set the number of iterations to 15 000 and the population size to 50 individuals). Furthermore, we also perform 30 runs to guarantee that the results do not depend on the stochasticity of the algorithm.

Numerical simulations and calibration of the model for the case of Italy

We present here the results obtained with our model for the case of Italy, as an example of application. The situation in that country, like in other European countries such as Spain, the United Kingdom and France, has been especially complicated since it has been strongly affected by the disease. The first two known cases of Coronavirus in Italy, a couple of Chinese tourists, were confirmed on January 30. The first known case of secondary transmission occurred in Codogno, Municipality of Lombardy, on February 18, 2020. After that, the pandemic grew and extended to other countries, which needed the establishment of strict control measures (see [52]).

In Section 4.1, we introduce the data and some model parameters specific to the case of Italy. Then, in Section 4.2, we show and analyze some results returned after the numerical experiments performed to illustrate the efficiency of our approach.

Considered data and model parameters

The information regarding the first cases in Italy is available in [52]. More precisely, the first two known cases of Coronavirus in Italy, a couple of Chinese tourists, were confirmed on January 30. On February 26, they were declared recovered. The first known case of secondary transmission occurred in Codogno, Municipality of Lombardy, on February 18, 2020. Thus, in the numerical simulations, we assume that the index case started his incubation period in Italy on =19 January 2020, as the earliest symptom onset of confirmed patients can be traced back to 30 January 2020 (reported on 31 January 2020) and we assume a delay of =10.5 days (rounded to 11 days) between the beginning of the incubation period and the symptoms onset. Thus, we set, , and all other compartments are set to 0.

The different control measures implemented in Italy are detailed in [52]. They started on February 2020, where a ’red zone’ was instituted on 11 municipalities with local transmission of SARS-CoV-2 infection, with lockdown to contain the emerging threat. Before this date we set and after this date we consider (to be estimated). On March 2020, lockdown was declared on the entire country. Thus, we assume with . On =22 March 2020 new additional measures regarding the containment and management of the epidemiological emergency from COVID-19, applicable on the whole national territory were adopted including the closure of non-essential or strategic production activities until May 3 (after several extensions). Consequently, we set , corresponding to the most restrictive social distancing control measures. Starting on May 2020, some production activities are authorized and we assume with . The Decree of 16 May 2020 establishes the end of regional travel restrictions and the reopening of production activities from May 2020. We set with . Finally, this Decree also governs the end of total travel restrictions (traveling to and from abroad can only be limited with measures derived from the legal system of the European Union and international obligations) and the reopening of production, commercial and social activities from =3 June to 31 July [53]. Hence, we consider with .

For Italy we do not include the compartment (i.e. people who die undetected pass from compartment to compartment ). This simplification is done since, according to [54], almost all undetected deaths affect to vulnerable persons that die quickly after developing clinical sign (e.g. aged persons in nursing homes). All the other compartments were used in this case study.

Furthermore, according to [55], noticeable changes in the delay of detection of infected persons (i.e. ) have not been observed during the course of the epidemic in Italy. Thus, we set .

The information by regions concerning the number of cases (i.e. people who enter each day in compartments or ) and COVID-19 deaths per day is daily reported in the website of the Italian Ministry of Health, (see [49]). Moreover, the Italian National Institute of Health supervises daily the data related to cases, deaths, tests carried out, etc. and the Italian civil protection publishes them completely here: In our simulations, we take , , , from the previous GitHub repository, from [56] and , from [48].

https://github.com/pcm-dpc/COVID-19

We use the multiobjective technique described above to identify the parameters taking and . We reduce the number of objective functions to two functions, and that combine the four time-series errors and , using some weights as follows: Thus, we solve the following two-objective problem: instead of solving a four objective one, that is more complex and time-consuming. We use (penalization to avoid the underestimation of the reported data in ) and (penalization to avoid the underestimation of the reported data at ).

Regarding the data denoted by , reported by the Italian authorities, after several numerical experiments, we decided to define as follows: Next, we solve the mono-objective problem described by (7) to identify the parameters .

The effective reproduction number given by (2) has been computed using the Composite Simpson’s rule.

Results

First, we use the methodology and techniques detailed in Sections 2, 3 to identify suitable model parameters for the case of Italy, validating the model proposed here. Then, we show how the model may be used to estimate the spread of the disease for different futures scenarios.

Validation of the model.

There is much uncertainty about but, as explained in Section 2.1 , we can first estimate and then for different values of (which is the instantaneous infection detected fatality ratio, iIdFR). That shows that, actually, the uncertainty of is related to the uncertainty of . Hence, we have considered different scenarios for the value of , using several estimations found in the literature for the Infection Fatality Ratio (IFR). We remind that the IFR is the expected value of the probability of dying for a person who is infected, i.e. the number of deaths due to the virus (detected or not) over the number of infections (detected or not). The IFR is one of the most critical features of the coronavirus disease, being quite difficult to estimate. More precisely:

According to [57], published on 4 May 2020, IFR=1.2%. Then, since on 5 May 2020 there were 29 315 detected deaths (see [58]) and 9154 estimated undetected deaths (see [21], [54]), an associated estimation of would be .

According to [59], published on 18 May 2020, IFR=1.91%. Then, an associated estimation of would be .

According to [60], published on 21 May 2020 but received on 4 April 2020, IFR=5.7%. Then, since on 4 April 2020 there were 15 362 detected deaths and on 14 April 2020 there were 6773 estimated undetected deaths (see [21], [22]), an associated estimation of would be .

We have used these three values of up to 21 July 2020.

The whole set of model parameters obtained when considering is reported in Table 1, Table 2. For the sake of simplicity, we only present here the parameters obtained with this particular value. Results for other values can be found at https://github.com/momat-ucm/T-SIR-T.

Table 1

Summary of all parameters used for COVID-19 in Italy when considering . A brief description (Description) of each parameter is given. Furthermore, if available, the range of the considered values (Value) and the reference in the literature (Ref.) are reported. The notation * means that this value is identified during the multiobjective or monoobjective optimization processes. The notation – means this parameter is omitted during this work. The notation & means that this parameter (depending on time) is estimated directly from the reported data. The notation means that the value of the parameter is set by the authors of this article, after some numerical experiments.

Notation	Value	Description	Ref.
N	60317000	Number of people in the country before starting the pandemic	[61]
γ_E	15.5	Transition rate of a person in compartment E (day−1)	[16]
γ_I(t)	15	Transition rate of a person in compartment I (day−1)at time t	[62]
γ_I_u(t)	19	Transition rate of a person in compartment I_u to compartment R_u (day−1) at time t	[62]
γ_H_R(t)	114.2729	Transition rate of a person in compartment H_R to compartment Q (day−1) at time t	*
γ_H_D(t)	15	Transition rate of a person in compartment H_D to compartment D (day−1) at time t	[62]
γ_I_D_u(t)	∞	Transition rate of a person in compartment I_D_u to compartment D_u (day−1) at time t	[54]
γ_Q,1	136.0450	Transition rate of a person in compartment Q to compartment R_d (day−1) at time t≤t_γ_Q	*
γ_Q,2	124.8646	Transition rate of a person in compartment Q to compartment R_d (day−1) at time t>t_γ_Q	*
t_γ_Q	2 May 2020	Date when the duration of a person in quarantine changes	*
d_g	0	Maximum number of days that d_I can be decreased	[55]
C_E	0.3806	Ratio between the disease contact rates β_E and β_I	[6]
C_u	0.3293	Ratio between β__I and the disease contact rate β_I	*
C_I_D_u	–	Ratio between the disease contact rates β_I_D_u and β_I	–
β_I,0	0.4992	Disease contact rate of a person in compartment I without control measures (day−1)	*
β_I_D_u,0	–	Disease contact rate of a person in compartment I_D_u without control measures (day−1)	–
β_H_R	&	Disease contact rate of a person in compartment H_R	&
β_H_D	&	Disease contact rate of a person in compartment H_D	&
t_0	19 Jan 2020	Initial time	[52]
ω_	1.4555	iIdFR when the implemented control measures are fully applied (in %)	[59]
δ_ω	0	to compute iIdFR when no control measures are applied ω¯=ω_+δ_ω	♢
ω_u,0	0.42	iIuFR when t=t_0 (in %)	*
ω_u,1	0.42	iIuFR when t=t_u,1 (in %)	♢
ω_u,2	0	iIuFR when t=t_u,2 (in %)	♢
t_u,1	3 May 2020	One day before the last report of d_r,u	[49]
t_u,2	4 May 2020	Date of the last report of d_r,u	[49]
θ(t)	&	Proportion of infected people that are detected by the authorities	&
p_0	0.7382	Ratio of the number of new detected infected people that will survive the disease and are hospitalized at time t_θ_0	[49]

Table 2

Summary of all parameters related to control measures used for COVID-19 in Italy when considering . A brief description (Description) of each parameter is given. Furthermore, if available, the range of the considered values (Value) and the reference in the literature (Ref.) are reported. The notation * means that this value is optimized during the multiobjective optimization process. The notation – means this parameter is omitted during this work.

Notation	Value	Description	Ref.
λ_0	19 Jan 2020	First day of application of the control strategy that was being used before t_0, if any	[52]
λ_1	23 Feb 2020	First day of application of the 1st control strategy (lockdown on 11 municipalities)	[52]
λ_2	11 Mar 2020	First day of application of the 2nd social distancing strategy (lockdown on the entire country)	[52]
λ_3	22 Mar 2020	First day of application of the 3rd social distancing strategy (closure of non-essential or strategic production activities)	[52]
λ_4	4 May 2020	First day of application of the 4th social distancing strategy (end of closure of non-essential or strategic production activities)	[52]
λ_5	18 May 2020	First day of application of the 5th social distancing strategy (end of regional travel restrictions and reopening of production activities)	[52]
λ_6	3 Jun 2020	First day of application of the 6th social distancing strategy (progressive end of total travel restrictions, reopening of production, commercial and social activities)	[53]
m_0	1	Intensity of the social distancing strategy used before t_0	[52]
m_1	1	Intensity of the control strategy that was being used in (t_0,λ_1]	[52]
m_2	0.5332	Intensity of the social distancing strategy used in (λ_1,λ_2]	*
m_3	0.1369	Intensity of the social distancing strategy used in (λ_2,λ_3]	*
m_4	0	Intensity of the social distancing strategy used in (λ_3,λ_4]	[6]
m_5	0.0549	Intensity of the social distancing strategy used in (λ_4,λ_5]	*
m_6	0.0577	Intensity of the social distancing strategy used in (λ_5,λ_6]	*
m_7	0.0578	Intensity of the social distancing strategy used in (λ_6,λ_7]	*
m^(t)	–	Intensity of the contact tracing control strategy	–
m^(s)	–	Intensity of the health system control strategy	–
κ_1	–	Efficiency of the social distancing strategy used in (λ_0,λ_1]	–
κ_2	62.0015	Efficiency of the social distancing strategy used in (λ_1,λ_2]	*
κ_3	0.5139	Efficiency of the social distancing strategy used in (λ_2,λ_3]	*
κ_4	0.0362	Efficiency of the social distancing strategy used in (λ_3,λ_4]	*
κ_5	90.6408	Efficiency of the social distancing strategy used in (λ_4,λ_5]	*
κ_6	14.6514	Efficiency of the social distancing strategy used in (λ_5,λ_6]	*
κ_7	52.0747	Efficiency of the social distancing strategy used in (λ_6,λ_7]	*
κ^(t)	–	Efficiency of the contact tracing control strategy	–
κ^(s)	–	Efficiency of the health system control strategy	–

In Figs. 2 and  3, we show the comparison of some of the outputs of the simulation run with the chosen parameters and the official real data, when . Similar figures have been obtained when considering and 3.9559% (obviously, with different parameters) and, thus, are not reported here.

Fig. 3

Comparison of some of the outputs of the simulation run with (Est.) and, when available, the Italian official reported data (Rep.). Top-left: New detected cases per day. Top-right: New detected deaths per day. Bottom-Left: Function modeling the social distancing measures. Bottom-Right: Effective reproduction number and the contribution of each infectious state (i.e. , , , and ) to this number. The vertical line corresponds to the first date when R.

More precisely, in Fig. 2, we show a very good fitting of the model outcomes with respect to four important time series reported by the Italian authorities regarding the pandemic: the cumulative number of cases and deaths, people in hospital, people in quarantine and recovered people. Furthermore, the model estimates that on 14 April 2020 and 5 May 2020 there were around 7509 and 9192 undetected deaths (persons in ), respectively. According to the Italian authorities (see [21], [54]), on 14 April 2020 and 5 May 2020 (the only two reported values) this number was estimated to be 6773 and 9154, respectively. This shows that our model is able to estimate well the biological behavior of the pandemic.

Next, in Fig. 3 (Top-Left) and (Top-Right), we show the evolution of the new detected cases and deaths, per day. When available, we compare those outputs with the real data. Again, we observe a good fitting between simulated results an observed data. This feature of the model is interesting as it allows to estimate the date and magnitude of the peaks of all those outcomes (also reported in this figure), which may help to plan some health system measures (such as increasing the number of beds in hospital).

In Fig. 3 (Bottom-Left) and (Bottom-Right), we show the evolution of the functions corresponding to the social distancing measures and the effective reproduction number, when . We observe that the effective reproduction number went below 1 after 14 March 2020, which should indicate the start of a decrease of the spread in Italy. This value is reached soon after the implementation of some intensive control measures (lockdown in the entire country) in Italy on 11 March 2020, showing that such restrictive control measures were necessary to control the pandemic (before this date the model estimates that R remains above 1). We also note that the initial value of the effective reproduction number, which is the basic reproduction number R, is equal to 8.0, which is a value higher than those generally reported in the literature (between 1.4 and 3.9, see e.g. [4], [63]). However, we also show in this figure the contribution of each infectious state (i.e. , , , and ) to the value of R (see Section 2.2 for more details). In particular, we observe that the compartment of undetected infectious people contributes to with a value of 3.9, which means that detected cases contribute to with a value of less than 4.1, which is consistent with the values usually found in the literature.

Next, in Fig. 4, we show the evolution of some of the outcomes of the model for different values of . First, in this figure (Top-Left) and (Top-Right), we present the evolution of R and the effect of the social distancing measures when considering , and 3.9559%. We observe that, before 11 March 2020, the magnitude and shape of those functions are affected by the value of . After this date, they are quite similar for the different values of . Focusing on the total number of cases (including undetected), presented in Fig. 4 (Bottom-Left), the model estimates that, taking into account the different values of considered previously, between 1 and 4 million people could have been infected by COVID-19 in Italy. Regarding the value of , depicted in Fig. 4 (Bottom-Right), we see that for the proportion of detection of new cases oscillates between 15% and 67%. For or , it is between 5.7% and 25%. From those results we see that, in all cases, the model estimates that a large amount of new cases have not been detected during this pandemic. However, the proportion of detection is highly dependent on the value of and large variations of the results are observed when this value changes. The results presented in Fig. 4, indicate the necessity to approximate well this critical parameter.

Fig. 4

Some of the outputs of the simulation run with and 3.9559% for the Italian case. Top-left: Effective reproduction number. Top-right: Function modeling the social distancing measures. Bottom-left: Cumulative number of total cases (including undetected cases). Bottom-right: .

Estimation of possible future scenarios

We considered just a few cases, to illustrate the possibilities that can be explored with the model to estimate the possible spread of the disease. We simulate the future scenarios described below until 1 September 2020. We recall that the meaning of parameters , and that appear below may be seen in Section 2.1. In all cases, we consider = 22 July 2020 (date of a potential change of control measures; taking into account that the date of our last data is 21 July 2020), . Since we cannot estimate for the future (it depends on the authority actions) by using the methodology described in Section 3.2, from 22 July 2020 on, we assume that remains equal to its value on 21 July 2020.

First, for , and , we consider . Some results of this simulation are reported in Fig. 5 (Top-Left). In the three cases we observe a slight increase in the number of cases.

Fig. 5

Outputs obtained when simulating possible future scenarios (using ). Top-left: New detected people when and .9145%, 1.4555% or 3.9559%, Top-right: Effective reproduction number when and , 0.25 or 0.4, Bottom-left: New detected cases per day, when and , 0.25 or 0.4, Bottom-right: New detected deaths per day, when and , 0.25 or 0.4.

Next, for , we simulate four future scenarios with: (i) ; (ii) ; (iii) , and (iv) . In Fig. 5 (Top-Right), (Bottom-Left) and (Bottom-Right), we show some outputs of the model for these four scenarios. As we can see, there is a remarkable possibility of new outbreaks in Italy if the control measures are relaxed too much. For instance, focusing on the effective reproduction number, it goes back over 1 when and , which would mean a second wave of uncontrolled outbreaks. In particular, if control measures are relaxed too much (e.g. ) the new cases and deaths may increase dramatically. This clearly indicates that, the pandemic is not over yet and efforts must continue, keeping enough control measures.

Conclusions and discussion

We have developed an epidemiological model for simulating the dynamics of COVID-19 in a considered territory. This model may be used by epidemiologists and policy makers for understanding the spread of the pandemic, study the efficiency of control measures and propose suitable ways to tackle the problem. It is intended to be as simple as possible, but with sufficient complexity to be able to simulate the biological and sociological mechanisms which influence the disease spread (including the three main types of control measures: social distancing, contact tracing and health system measures), to fit real data and to simulate possible future scenarios. Of course, for particular cases the model may be simplified (for instance removing the compartment, as we showed for the case of Italy, or the compartment, if the territory does not keep asymptomatic people in quarantine at home).

The model can serve as a basis for studies including higher complexity. Indeed, it can be modified to introduce structures of age, sex, etc., stochasticity, levels of uncertainty, mobility between territories, spatial structure, etc.

In the model that we propose, there are several important points to highlight that, to the best of our knowledge, have not been previously considered:

We have included a compartment for undetected deaths mainly happening in nursing homes or residences for disabled people. The number of these deaths and the associated infections does not seem to be negligible, which makes necessary to take them into account in the model.

Some dynamic coefficients are based on reported time series, which allows the model to reproduce faithfully the spread of the disease, fitting very well the reported data.

The model includes the effects of the different control measures described above.

The model takes into account quarantine and hospitalized people in order to estimate the number of hospital beds that are needed.

We point out that the model, in this current simple form, has the following limitations:

Epidemics/pandemics are phenomena that respond to complex processes and all epidemiological models are simplifications of this complex reality. Furthermore, the model is based on the interactions of human beings, which do not behave in reasonably predictable ways like molecules, cells or particles (see [35]).

The prediction of its evolution is very complicated because of the characteristics of the system of differential equations governing the model. In particular, in the exponential growth phases, small variations in some parameters can lead to large variations in model outcomes (see, [6]). Moreover, there are data quality problems: a) uncertainty about the characteristics of SARS-CoV-2 (new virus); (b) changes in criteria to count cases; (c) uncertainty about undetected cases (and, therefore, the mortality rate), etc.

Since there is no clear scientific evidence of the effect of humidity and temperature (see, e.g. [31], [32]), we have not included these two factors in our model (if such evidence is demonstrated, these variables should be included in the model).

It is assumed that the spatial distribution of the population is homogeneous in the territory under study (if this is not a suitable assumption, the territory can be divided into smaller ones).

In this paper we have only considered the spread of disease within a territory, where only local transmission and known imported/exported cases are modeled. The between-country spread has not been modeled in this work, but can be considered in the future employing ideas from the BeCoDiS model (see [15]).

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.

Table 3

Summary of the parameters used to model the control measures corresponding to the days after the 21 July 2020. Time intervals with different measures are indexed with .

Date of application	Efficiency	Intensity	Type	Detail
(λ_i−1)	(κ_i)	(m_i)
1 Aug 2020	62.0015	0.3	Relaxation	“New normality”
20 Sep 2020	62.0015	0.53	Relaxation	Back from vacation,
				school, work. Regional
				elections
3 Nov 2020	62.0015	0.3	Restriction	Curfew
13 Nov 2020	62.0015	0.22	Restriction	More regions with
				restrictive measures
23 Dec 2020	62.0015	Different	Different	Christmas
		scenarios	scenarios
7 Jan 2020	62.0015	Different	Different	End of Christmas
		scenarios	scenarios