key: cord-0784323-gsxewkw4 authors: Ramos, A.M.; Ferrández, M.R.; Vela-Pérez, M.; Kubik, A.B.; Ivorra, B. title: A simple but complex enough [Formula: see text]-SIR type model to be used with COVID-19 real data. Application to the case of Italy date: 2021-01-01 journal: Physica D DOI: 10.1016/j.physd.2020.132839 sha: 660536713221bee9446123f7035617854bde8319 doc_id: 784323 cord_uid: gsxewkw4 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 type 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 [Formula: see text]-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. The outbreak of novel coronavirus disease 2019 continues to spread rapidly around the world [60] . 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 [63] when the disease was confirmed in more than 118,000 cases reported globally in 114 countries. Although some authors (see, e.g., [37] ) 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 China and Senegal all infected people are treated in hospitals (or other kind of health facilities) and, therefore, compartment Q is not necessary [33] . 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. 2. 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. 3. 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. 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. According to the known characteristics of the COVID-19 pandemic, we assume that each person is in one of the following compartments (see [2, 47, 61, 57] ): • Susceptible (denoted by S): The person is not infected by the disease pathogen. • Exposed (denoted by E): 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., [2] ), 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 I. • Infectious (denoted by I): 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., [61] ). 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 other compartment, I u or I Du ). • Infectious but undetected that will survive the disease (denoted by I u ): After being in the compartment I, 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 R u . • Infectious but undetected that will die (denoted by I Du ): After being in the compartment I, 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 D u . We remark that vulnerable people died outside hospitals during the first wave of the pandemic in some countries (such as Italy or Spain; see [32, 30, 24, 14, 41] ). 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 [32, 30] ), 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 [41] ). 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 I Du . We note that the results reported later in Section 4 indicate that this choice is adequate. • Quarantine (denoted by Q): The person is in quarantine at home or some special facilities (different from hospitals) during a fixed period. For example, according to [9, 29] , 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 I) and the model takes that into account. Furthermore, according to [19] , 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 Q 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 R d . 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 H R ): The person is in hospital and (according to [19, 51, 58] ) can still infect other people (in particular healthcare workers). At the end of this state, a person passes to the compartment Q (if implemented in the territory under study; otherwise the person passes directly to the compartment R d ). • Hospitalized that will die (denoted by H D ): The person is hospitalized and can still infect other people. At the end of this state, a person passes to the compartment D. • Dead by COVID-19 (denoted by D) after being at hospital. • Dead by COVID-19 but undetected (denoted by D u ): 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 R d ): 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 [47] , 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 R u ): 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 H R and H D 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 Figure 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). The authorities may apply various control measures in order to control the COVID-19 spread (see [6] ): • 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 [64] ). 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 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 I. 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 t 0 ≥ 0, considering that only susceptible people live in the territories that are free of the disease, whereas the number of people in compartments S, E, I, I u ,I Du , H R , H D , R d , R u , D and D u of the 6 J o u r n a l P r e -p r o o f infected territories are set to their corresponding values. Then, for the time interval [t 0 , t 0 + T max ], with T max ∈ IN being the maximum number of simulation days, the model is applied. If at the end of a simulation day t there is less than 1 person in each compartment E, I, I u , I Du , H R and H D in all the territory, the simulation is stopped. Else, the simulation is stopped when t = t 0 + T max . 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 Figure 1 . We assume that people in a territory are characterized to be in one of the compartments described above: S, E, I, 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 t 0 , 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 [35] with the Be-CoDiS model. There are some studies about the relevance of humidity and temperature for the spread of COVID-19. In [40] , 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 [56] 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 J o u r n a l P r e -p r o o f Journal Pre-proof of ordinary differential equations: (1) In this system we have used, as stated above, that people in compartment Q cannot infect other people. To include the possibility of infection, we would only need to add the term β Q Q inside the parenthesis of the first two equations. In (1): • N is the number of people in the territory before the start of the pandemic. • ω(t) ∈ [0, 1] is the instantaneous infection detected fatality ratio (iIdFR) in the territory, at time t: the proportion of new detected people (after being in compartment I) that will die of COVID-19 after 1 γH D (t) days, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time. More precisely, where c m (t) (respectively, c m,u (t)) is the model cumulative number of COVID-19 detected (respectively, undetected) cases at day t (see more details about them below, when describing the outputs of the model). • ω u (t) ∈ [0, 1] is the instantaneous infection undetected fatality ratio (iIuFR) in the territory, at time t: the proportion of new undetected people (after being in compartment I) that will die of COVID-19 after We point out that ω(t)+ω u (t) ∈ [0, 1] is the instantaneous infection fatality ratio (iIFR), at time t: 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 above (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, ω + ω u . 1] is the proportion of the number of new infectious people that are detected and documented by the authorities (after being in compartment I), per unit time, compared to the total number of new infectious people (detected or undetected), per unit time, at time t. More precisely, . We point out that is the instantaneous case fatality ratio (iCFR), at time t: the proportion of the number of new detected people (after being in compartment I) that will die of COVID-19 after 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, ω θ . • β E , β I , β Q , β Iu , β IDu , β HR , β HD ∈ [0, ∞) are the disease contact rates (day −1 ) of a person in the respective compartments. They may change with time, due to the control measures. • γ I (t) ∈ (0, +∞) is the transition rate (day −1 ) from compartment I to compartments I u , I Du , Q, H R or H D , at time t. Similarly, γ E , γ Iu (t), γ ID u (t), γ Q (t), γ HR (t), and γ HD (t) ∈ (0, +∞) denote the transition rates (day −1 ) from compartments E, I u , I Du , Q, H R and H D to compartments I, R u , D u , R d , Q and D, respectively, at time t. For some particular territories, some compartments could be excluded from the model. • p is the ratio of the number of new detected infected people that will survive the disease and are hospitalized (they go to compartment H R ) over the total number of new detected infected people that will survive the disease (and may also do quarantine, compartment Q, without going to hospital). • τ 1 (t) and τ 2 (t) 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 [35] . For the sake of simplicity we have included those terms only in the second equation, corresponding to compartment E, 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. . These values must be non-negative and It is standard (see, e.g., [34] , [36] ) to prove that S(t), E(t), We remark that system (1) is an epidemic-type model (as described, e.g., in [22] ), 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 S + R d + R u + D u + D = N 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: Q 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 Q, then p = 1, For the numerical simulations presented in this paper, the first eight equations of system (1) were numerically solved with the classic 4 th 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). Here, we present some of the outputs that can be obtained from the mathematical model, used to analyze the results of the simulations: • c m (t): The model cumulative number of COVID-19 detected cases at day t, which is given by and can be also computed as follows: • c m,n (t): The number of new cases per day at day t, which is given by • h m (t): The model cumulative number of COVID-19 detected infected healthcare workers at day t, which is given by • c m,u (t): The model cumulative number of COVID-19 undetected cases at day t, which is given by Journal Pre-proof and can be also computed as follows: The model cumulative number of detected deaths (due to , at day t, which is given by D(t). • d m,u (t): The model cumulative number of undetected deaths, at day t, which is given by D u (t). • R e : 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 [5] . 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 t =0 no control measures have been applied, R e (0) = R 0 , where R 0 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 e (t) < 1. 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 t, 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 0 can be also considered. For instance, one can use the so-called "Next generation Matrix" method to approximate its value (see [54] ). However, as pointed in [7] , this approximation does not correspond accurately to the biological definition of R 0 (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 e should be used with caution by policymakers, without relying only on that value (see [1] ). Note that we can estimate the particular contribution of each infectious compartment (i.e. E, I, I u , I Du , H R and H D ) to the value of R e and R 0 by computing only their respective integrals in (2). • Hos(t): 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. • Hos n (t): The number of new hospitalized people per day at day t, which is given by • MHos: The maximum number of hospitalized people at the same time (in territory i) during the time interval [t 0 , T ]. It is computed as: Hos(t). This number can help to estimate and plan the number of clinical beds needed to treat all the COVID-19 cases. • Γ E (t) and Γ Iu (t): The number of people infected during the time interval [t 0 , T ], by contact with people in compartment E and I u , 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. 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. 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 β E,0 , β I,0 , β Iu,0 , β ID u ,0 , β HR,0 , β HD,0 are the respective disease contact rates when no control measures are taken and m E , m I , m Iu , m ID u , m HR , m HD are functions simulating de control measures, as described below. Regarding β Iu,0 , we may consider nonincreasing , if θ ∈ (ω, 1), where ω is a fatality rate value that is explained below and β I and β I are suitable lower and upper bounds, respectively. For the sake of simplicity, here we take β I,0 = β I,0 . Furthermore, following the idea proposed in [33, 35, 39] , we assume that there exists a relationship between the contact rates β I,0 , β E,0 , β ID u ,0 and β I,0 . More precisely, following [33] , we consider that people in compartments E, I u and I Du are less infectious than people in compartment I (due to their probably lower virus load or isolation measures, see [2, 39, 61] ). Thus, we consider that β E,0 = C E β I,0 and β I,0 = C u β I,0 and β ID u ,0 = C ID u ,0 β I,0 , where C E , C u , C ID u ∈ [0, 1]. A similar approach can be followed for β HR and β HD assuming, for instance, that β HR,0 = β HD,0 = C H β I,0 . Alternatively (as we do here for the case of Italy) we can compute good estimates for β HR and β HD,0 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 m E (t), m I (t), m Iu (t) ∈ [0, 1] (%) represent the reduction of some of the disease contact rates, at time t, due to the efficiency of the control measures corresponding to social distancing applied to the compartments E, I and I u , respectively (see, e.g., [46] ). 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 [33] ) We may also define in a similar way functions m (t) and m (s) related with contact tracing and health system control measures, respectively, which include similar parameters denoted by m j 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 c i instead of m i , as follows: • If we want m i+1 ≤ m i , then we look for a suitable parameter c i+1 ∈ [0, 1] such that m i+1 = c i+1 m i . • If we want m i+1 ≥ m i , then we look for a suitable parameter c i+1 ∈ [0, 1] such that m i+1 = m i + c i+1 (1 − m i ). A similar approach can be followed by using new parameters c Function θ is related with the contact tracing measures, since the higher those measures are at time t, the closer θ(t) is to 1. The other function involved with the contact tracing measures is m (t) , which is used below for estimating the transition rates. The function m (s) is used to compute the fatality ratios and to define m HR (t) = m HD (t) = m (s) (t) ∈ [0, 1], representing the reduction of the disease contact rates β HR,0 , β HD,0 , at time t, 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 β HR and β HD using historical data (which will be shown below). . 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., [3, 53] ). Additionally, we assume that , if t = t u,n constant after the last value t u,n linear otherwise, where n ∈ N ∪ {0} and the corresponding values t i , ω u,i are parameters to be chosen for each territory, depending on its particular circumstances. -Transition rates: • According to [39] , the transition rate γ E from E to I depends only on the disease and, therefore, is considered constant. • Some other transition rates could depend on time. For instance, according to [39, 38] , depending on the studied territory, the value of γ I (t) 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 γ Iu (t), γ ID u (t), γ Q (t), γ HR (t) and γ HD (t) can be decreased (e.g. people with symptoms stay under observation during more time). In order to parametrize these functions, we denote by d I , d Iu , d ID u , d Q , d HR and d HD the initial duration (given as a suitable average) in days of a person in compartment I, I u , I Du , Q, H R and H D , respectively, without the application of special contact tracing measures. Then, following [35] , we may use the following functions: where g(t) = d g (1 − m (t) (t)) represents the decrease of the duration of d I due to the application of contact tracing measures, at time t; and d g is the maximum number of days that d I 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 d g = 0. -Ratio p: We describe it by using the following function depending on the value of θ − ω (we remark that other alternatives are also possible): Here p 0 = p(t θ0 ) and ω 0 = ω(t θ0 ), where t θ0 and θ 0 are described below (when we show how to estimate θ(t)). Indeed, we assume that, if θ − ω > θ 0 − ω 0 all new detected infected people, when comparing to the case of θ 0 and ω 0 , will go to state Q, as they should present light symptoms. Furthermore, we use a linear interpolation between the case θ − ω = 0 (with p = 1, i.e., only persons with severe symptoms are detected) and θ − ω = θ 0 − ω 0 (with p = p 0 ). 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 θ(t), has not been parameterized above. We will show below how to deal with it. 14 J o u r n a l P r e -p r o o f Journal Pre-proof 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 t. • d r (t): The reported cumulative number of COVID-19 deaths, at time t. • h r (t): The reported cumulative number of COVID-19 infected healthcare workers, at time t. • τ 1,r (t): The reported number of known imported cases on day t. 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). • τ 2,r (t): The reported number of known exported cases on day t. 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). • d r,u (t): The reported cumulative number of estimated undetected COVID-19 deaths, at time t. • Hos r (t): The reported number of people in hospital (i.e. Hos r ), at time t. • q r (t): 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., [17] ). Now, with that information we may estimate some of the model parameters and functions, as follows: -Ratio θ(t) of new infected people that are detected: To estimate θ(t), first, we compute an estimation of iCFR ( ω(t) θ(t) ) as follows: , r is the minimum natural number such that c r (t) − c r (t − r) = 0 and the previous quotient is > 1, where t iCFR ≥ 7 days is the first date for which c r (t) − c r (t − 1) = 0 and d r (t + 1 γH D ) = 0. 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., [23] ). 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 ω CFR (t): To avoid unrealistic situations (due to poor data), if ω CFR (t) ≤ 0.01 (which is assumed too low, see, e.g., [53] ) for some t > t iCFR , then we extend the average of the last week used to compute this ω CFR (t) by including the previous days until getting ω CFR (t) > 0.01. Finally, we define θ(t) by: and t θ0 = t iCFR + 6 days to guarantee that ω CFR (t θ0 ) 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 ω(t) θ(t) . In this case, we set θ(t) 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 θ(t). -Contact rates β HR (t) and β HD (t): We use the following estimation of the fraction of infected healthcare workers (assumed to be due to contact with hospitalized people, see, e.g., [58] ): To compute η, we first define η as and this quotient is ≤ 1, and t η is the first date for which h r t η + 1 γE + 1 γI = 0 is reported. We point out that, in the previous formula, we have made a distinction about the knowledge of the values of h r , because they are rarely given every day. We remark that we are taking into account the delay of 1 γE + 1 γI days between the infection of a person and its detection by authorities. Thus, we set η(t) to a constant value for the last 1 γE + 1 γI 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 β HR = β HD , we have that In the previous expression, the coefficients depend on time t. -Other unknown parameters: We can identify them with a multiobjective technique, as described below. 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 S ⊂ R n the feasible region given by the cartesian product of all those intervals, where n is the number of unknown parameters. Once S is defined, we follow the methodology proposed in [15] . 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 Ω ∈ S. Then, given two weight coefficients ρ 1 > 0 and ρ 2 > 0 and a temporal series X, we define the function that measures the error between the reported data, which is denoted by X r , and the model output denoted by X Ω m . Notice that, as done in previous works [15, 33] , 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 [50] , we use the L 2 norm in [t 0 , t 0 + T max ], which is The number of functions f X (Ω) 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. c r and d r ), see, e.g., [62] . In some cases, they also include other data as, for instance, the number of people in hospital (i.e. Hos r ), in quarantine at home (i.e. q r ), or the cumulative number of undetected deaths (i.e. d r,u ) (see, e.g., [44] ). Then, we formulate a multiobjective problem that may have as many objective functions as available temporal series. For instance, we could solve where S is the feasible region (i.e. upper and lower bounds of each parameter) and n 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 [16, 17, 52] . 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 [17] . 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 [52] ) 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 [52] 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 15000, 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 [52] . 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 equations of system 1, which, as said previously, are not coupled with the other equations, there is a parameter, namely γ Q , 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. q r ), the parameter γ Q can be estimated separately as a final step of the identification methodology. To do that, we assume that γ Q can be parameterized by n ∈ N parameters denoted by Φ = {Φ 1 , · · · , Φ n } ∈ S ⊂ R n (see next section for a particular example). Then, we formulate the following mono-objective problem: where f q (Φ) measures the error (as defined by (5)) between the reported data q r and q Φ m , the number of people in quarantine returned by the model when considering γ Q 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 15000 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. 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 [42] ). 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. The information regarding the first cases in Italy is available in [42] . 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 t 0 =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 d I + d E =10.5 days (rounded to 11 days) between the beginning of the incubation period and the symptoms onset. Thus, we set, E(t 0 ) = 1, S(t 0 ) = N − 1 and all other compartments are set to 0. The different control measures implemented in Italy are detailed in [42] . They started on λ 1 =23 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 m 0 = m 1 = 1 and after this date we consider m 2 ∈ [0, 1] (to be estimated). On λ 2 =11 March 2020, lockdown was declared on the entire country. Thus, we assume m 3 = c 3 × m 2 with c 3 ∈ [0, 1]. On λ 3 =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 m 4 = 0, corresponding to the most restrictive social distancing control measures. Starting on λ 4 = 4 May 2020, some production activities are authorized and we assume m 5 = c 5 × m 2 with c 5 ∈ [0, 1]. The Decree of 16 May 2020 establishes the end of regional travel restrictions and the reopening of production activities from λ 5 =18 May 2020. We set m 6 = m 5 + c 6 × (1 − m 5 ) with c 6 ∈ [0, 1]. 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 λ 6 =3 June to 31 July ( [43] ). Hence, we consider For Italy we do not include the compartment I Du (i.e. people who die undetected pass from compartment I to compartment D u ). This simplification is done since, according to [31] , 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 [26] , noticeable changes in the delay of detection of infected persons (i.e. d I ) have not been observed during the course of the epidemic in Italy. Thus, we set d g = 0. The information by regions concerning the number of cases (i.e. people who enter each day in compartments H R or H D ) and COVID-19 deaths per day is daily reported in the website of the Italian Ministry of Health, (see [44] ). 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: https://github.com/pcm-dpc/COVID-19 In our simulations, we take c r (t), d r (t), Hos r (t), q r (t) from the previous github repository, h r (t) from [28] and τ 1,r (t), τ 2,r (t) from [62] . We use the multiobjective technique described above to identify the parameters Ω = (β I,0 , κ 2 , κ 3 , κ 4 , κ 5 , κ 6 , κ 7 , c 2 , c 3 , c 5 , c 6 , c 7 , p 0 , C u , ω u,0 , 1/γ HR ) ∈ R 16 , taking ρ 1 = 10 and ρ 2 = 100. We reduce the number of objective functions to two functions, f 1 (Ω) and f 2 (Ω) that combine the four time-series errors f c (Ω), f d (Ω), f Hos (Ω) and f du (Ω), using some weights as follows: and Thus, we solve the following two-objective problem: instead of solving a four objective one, that is more complex and time-consuming. We use ρ 1 = 10 (penalization to avoid the underestimation of the reported data in [t 0 , t 0 + T max )) and ρ 2 = 100 (penalization to avoid the underestimation of the reported data at t 0 + T max ). Regarding the data denoted by q r (t), reported by the Italian authorities, after several numerical experiments, we decided to define γ Q (t) as follows: Next, we solve the mono-objective problem described by (7) to identify the parameters Φ = (γ Q,1 , γ Q,2 , t γQ ) ∈ [0, 1] × [0, 1] × [0, +∞). The effective reproduction number given by (2) has been computed using the Composite Simpson's rule. First, we use the methodology and techniques detailed in Sections 2 and 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. There is much uncertainty about θ but, as explained in the Section 2.1 , we can first estimate ω(t) θ(t) and then θ(t) 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 [55] , published on 4 May 2020, IFR=1.2%. Then, since on 5 May 2020 there were 29315 detected deaths (see [13] ) and 9154 estimated undetected deaths (see [31, 32] ), an associated estimation of ω would be ω = 1.2% 29315 29315+9154 = 0.9145%. • According to [21] , published on 18 May 2020, IFR=1.91%. Then, an associated estimation of ω would be ω = 1.91% 29315 29315+9154 = 1.4555%. • According to [11] , published on 21 May 2020 but received on 4 April 2020, IFR=5.7%. Then, since on 4 April 2020 there were 15362 detected deaths and on 14 April 2020 there were 6773 estimated undetected deaths (see [30, 32] ), an associated estimation of ω would be ω = 5.7% 15362 15362+6773 = 3.9559%. J o u r n a l P r e -p r o o f We have used these three values of ω up to 21 July 2020. The whole set of model parameters obtained when considering ω = 1.4555% is reported in Tables 1 and 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. In Figures 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 ω = 1.4555%. Similar figures have been obtained when considering ω = 0.9145% and 3.9559% (obviously, with different parameters) and, thus, are not reported here. More precisely, in Figure 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 D u ), respectively. According to the Italian authorities (see [32, 31] ), 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 Figure 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 Figure 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 ω = 1.4555%. 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 e remains above 1). We also note that the initial value of the effective reproduction number, which is the basic reproduction number R 0 , 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. [8, 39] ). However, we also show in this figure the contribution of each infectious state (i.e. E, I, I u , H R and H D ) to the value of R e (see Section 2.2 for more details). In particular, we observe that the compartment I u of undetected infectious people contributes to R 0 with a value of 3.9, which means that detected cases contribute to R 0 with a value of less than 4.1, which is consistent with the values usually found in the literature. Next, in Figure 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 e and the effect of the social distancing measures when considering ω = 0.9145%, ω = 1.4555% 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 Figure 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 Figure 4 (Bottom-Right), we see that for ω = 3.9559% the proportion of detection of new cases oscillates between 15% and 67%. For ω = 0.9145% or ω = 1.4555%, 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 Figure 4 , indicates the necessity to approximate well this critical parameter. 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 λ i , κ i and m i that appear below may be seen in Section 2.1. In all cases, we consider λ 7 = 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), κ 8 = κ 2 = 62.0015. Since we cannot estimate θ for the future (it depends on the authority actions) J o u r n a l P r e -p r o o f Journal Pre-proof Table 1 : Summary of all parameters used for COVID-19 in Italy when considering ω = 1.4555%. 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 Transition rate of a person in compartment E (day −1 ) [64] γ I (t) 1 5 Transition rate of a person in compartment I (day −1 ) at time t [25] γ Iu (t) 1 9 Transition rate of a person in compartment I u to compartment R u (day −1 ) at time t [25] γ HR (t) 1 Transition rate of a person in compartment H R to compartment Q (day −1 ) at time t * γ HD (t) 1 5 Transition rate of a person in compartment H D to compartment D (day −1 ) at time t [21] δ ω 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 [44] -Efficiency of the health system control strategy 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 ω = 0.9145%, ω = 1.4555% and ω = 3.9559%, we consider m 8 = 0.25. Some results of this simulation are reported in Figure 5 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 type 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 I Du compartment, as we showed for the case of Italy, or the Q 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 In the version of this article initially submitted to Physica D: Nonlinear Phenomena, the case of Italy studied in Section 4 was written using official data up to 21 July 2020 to validate the model and analyse several scenarios up to 1 September 2020. In the current version, written later as a revised version of the submitted paper, we have kept everything that was computed using those data (including the values of the parameters). Furthermore, following the reviewers suggestions, in this Annex we show an update of the situation in Italy, using official data up 21 December 2020, to compare with our simulations. It is important to remark that in all the simulations shown in this Annex, we keep the same parameter values as those given in the version of the paper initially submitted to the journal (with data up to 21 July 2020, as mentioned above). In all cases, we consider ω = 1.4555% and we estimate θ by using the methodology described in Section 3.2 until 21 December 2020 assuming that, since then, θ remains equal to its value on this day. Between 14 and 15 August 2020, the cumulative number of deaths reported by the Italian authorities shows a sudden increase, because they updated unreported deaths that occurred in previous weeks. According to their weekly reports (see [27] ), from 10 to 23 August there were 14 new deaths. Thus, 218 deaths must be smoothly distributed before the 10 August as done in [33] for the case of China. First, we distribute the 14 deaths from 10 to 23 August and, then, we assign the 218 deaths between 22 July and 23 August, avoiding the modification of the data that we had already used in the paper until 21 July. Besides the values of the model parameters shown in Section 4, in Table 3 we show the values of the parameters that we have used to model the control measures (as implemented for the Italian Government) corresponding to the days after the 21 July 2020. In Figure 6 we show some results corresponding to five different scenarios up to 31 January 2021: • Red scenario: More restrictive control measures are implemented from 23 December to 31 January. • Yellow scenario: More restrictive control measures are implemented from 23 December to 6 January. From 7 January, the measures come back to the ones before 23 December. • Purple scenario: The control measures are maintained without changes. • Green scenario: The control measures are relaxed for Christmas from 23 December to 6 January. From 7 January, more restrictive control measures are implemented. • Cyan scenario: The control measures are relaxed for Christmas from 23 December. These measures continue during January, i.e., no restrictions are implemented after Christmas. As we can see, Italy could be on 22 December at a crucial moment where any changes in control measures could have notable consequences. The model estimates that on 22 December the effective reproduction number R e could be slightly higher than 1 and the intensity of the control measures could be of 0.22 (considering 1 the absence of measures before the pandemic). If they tighten the measures from 23 December, as the Italian government has announced, the effective reproduction number could drop below 1 and the pandemic would be controlled (see the red scenario in Figure 6 ). However, if the measures are relaxed too soon in January, the cases and deaths would increase again before the end of January (see the yellow scenario in Figure 6 ). As we can see in the purple scenario, if no changes are implemented in the measures from 22 December to 31 January, the new cases and deaths may increase. Finally, if the measures are relaxed on 23 December, there would be a third wave of uncontrolled outbreaks whose magnitude depends on the restrictive measures they apply after Christmas (see the green and cyan scenarios in Figure 6 ). A guide to R -the pandemic's misunderstood metric Presymptomatic SARS-CoV-2 infections and transmission in a skilled nursing facility 12 =0.15, m 13 =0 12 =0.15, m 13 =0 12 =0.35, m 13 =0 Outputs obtained when simulating possible future scenarios when ω = 1.4555%. Top-left: Function modeling the social distancing measures Estimating the infection fatality rate among symptomatic COVID-19 cases in the United States The geography of COVID-19 spread in Italy and implications for the relaxation of confinement measures COVID-19 control in China during mass population movements at New Year. The Lancet The failure of R 0 Assessment of the sars-cov-2 basic reproduction number, r0, based on the early phase of covid-19 outbreak in italy Centers for Disease Control and Prevention. Transmission of SARS-COV-2 Infections in Households -Tennessee and Wisconsin. Mobidity and Mortality Weekly Report A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic The COVID-19 infection in italy: A statistical study of an abnormally severe disease Understanding Infectious Disease Dynamics Aggiornamento casi COVID-19 European Centre for Disease Prevention and Control. Surveillance of COVID-19 at long-term care facilities in the EU/EEA A multi-objective approach to estimate parameters of compartmental epidemiological models. Application to Ebola Virus Disease epidemics Highperformance computing for the optimization of high-pressure thermal treatments in food industry Preference-based multiobjectivization applied to decision support for High-Pressure Thermal processes in food treatment Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy European and US lockdowns and second waves during the COVID-19 pandemic Estimating the global infection fatality rate of COVID-19. medRxiv The basic epidemiology models: models, expressions for R 0 , parameter estimation, and applications Global Stability for Delay SIR and SEIR Epidemic Models with Nonlinear Incidence Rate Vigilancia de los excesos de mortalidad por todas las causas Characteristics of SARS-CoV-2 patients dying in Italy Report based on available data on Epidemia COVID-19, aggiornamento nazionale Epidemia COVID-19, aggiornamento nazionale Recommendation for people in self-isolation and for the family members assisting them Survey nazionale sul contagio COVID-19 nelle strutture residenziali e sociosanitarie Survey nazionale sul contagio COVID-19 nelle strutture residenziali e sociosanitarie Impatto dell'epidemia COVID-19 sulla mortalità totale della popolazione residente primo quatrimestre 2020 Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. the case of China Stability and sensitivity analysis of Be-CoDiS, an epidemiological model to predict the spread of human diseases between countries. Validation with data from the 2014-16 West African Ebola Virus Disease epidemic Be-CoDiS: A mathematical model to predict the risk of human diseases spread between countries. Validation and application to the 2014 Ebola Virus Disease epidemic Parameter estimation for a mathematical model predicting the COVID-19 spread in the Àrea Metropolitana de Barcelona Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet Infectious Diseases Transmission dynamics of 2019 novel coronavirus (2019-ncov). bioRxiv, 2020 Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2) The role of absolute humidity on transmission rates of the COVID-19 outbreak Informe del GTM sobre el impacto de la COVID-19 en las personas mayores, con especial énfasis en las que viven en residencias Covid-19 -Situazione in Italia Fase 2, nuovo decreto sulle riaperture: ecco cosa prevede COVID-19: Development of a robust mathematical model and simulation package with consideration for ageing population and time delay for control action and resusceptibility Mathematical assessment of the impact of nonpharmaceutical interventions on curtailing the 2019 novel Coronavirus Detection of sars-cov-2-specific humoral and cellular immunity in COVID-19 convalescent individuals A complete analysis of the epidemiological scenario around a SARS-CoV-2 reinfection: previous infection events and subsequent transmission. Preprint (Version 1) available at Research Square COVID-19 and Italy: what next? Health Policy Machine Learning with R, the tidyverse, and mlr Short-term Forecasts of the COVID-19 Epidemic in Guangdong and Zhejiang A preference-based evolutionary algorithm for multiobjective optimization: the weighting achievement scalarizing function genetic algorithm Estimating the infection and case fatality ratio for COVID-19 using age-adjusted data from the outbreak on the Diamond Princess cruise ship. medRxiv Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Using mobility to estimate the transmission intensity of COVID-19 in italy: a subnational analysis with future scenarios Temperature Significantly Change COVID-19 transmission in 429 cities Unique epidemiological and clinical features of the emerging 2019 novel coronavirus pneumonia (COVID-19) implicate special control measures Reasons for healthcare workers becoming infected with novel coronavirus disease 2019 (COVID-19) in China World Bank Open Data World Health Organization. Coronavirus disease (COVID-2019) Situation reports World Health Organization. Coronavirus disease (COVID-2019) Situation reports 73 World Health Organization. Coronavirus disease (COVID-2019) Situation reports Director-general's opening remarks at the media briefing on COVID-19 -11 Report of the WHO-China Joint Mission on Coronavirus Disease This research was partially supported by the Spanish Government under projects PID2019-106337GB-I00, MTM2015-64865-P and MTM2017-89423-P; the research group MOMAT (Ref. 910480) supported by the Complutense University of Madrid; and the Regional Government of Andalusia under project P12-TIC301, partially financed by the European Regional Development Fund (ERDF). Also thanks to M. Ramos-Rubio for his help with the manuscript. Est. Re =1.4555% Est. Re E =1.4555%Est. Re I =1.4555%Est. Re Iu =1.4555%Est. Re HR =1.4555%Est. Re HD =1.4555% 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:1. 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. 4 . 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:1. 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 [22] ).2. 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, [33] ). 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.3. Since there is no clear scientific evidence of the effect of humidity and temperature (see, e.g. [40, 56] ), we have not included these two factors in our model (if such evidence is demonstrated, these variables should be included in the model).4. 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).5. 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 [35] ).