key: cord-0754428-wh4ald1r authors: Carcione, Jose' M title: A simulation of a COVID-19 epidemic based on a deterministic SEIR model date: 2020-04-24 journal: nan DOI: 10.1101/2020.04.20.20072272 sha: 8c0613318fde14d3522efae9f3e7b0dd063098c8 doc_id: 754428 cord_uid: wh4ald1r An epidemic disease caused by a new coronavirus has spread in Northern Italy with a strong contagion rate. We implement an SEIR model to compute the infected population and number of casualties of this epidemic. The example may ideally regard the situation in the Italian Region of Lombardy, where the epidemic started on February 25, but by no means attempts to perform a rigorous case study in view of the lack of suitable data and uncertainty of the different parameters, namely, the variation of the degree of home isolation and social distancing as a function of time, the number of initially exposed individuals and infected people, the incubation and infection periods and the fatality rate. First, we perform an analysis of the results of the model, by varying the parameters and initial conditions (in order the epidemic to start, there should be at least one exposed or one infected human). Then, we consider the Lombardy case and calibrate the model with the number of dead individuals to date (April 19, 2020) and constraint the parameters on the basis of values reported in the literature. The peak occurs at day 37 (April 1) approximately, when there is a rapid decrease, with a reproduction ratio R0 = 3 initially, 1.38 at day 22 and 0.64 after day 35, indicating different degrees of lockdown. The predicted death toll is almost 14000 casualties, with 2.4 million infected individuals at the end of the epidemic. The incubation period providing a better fit of the dead individuals is 4.25 days and the infection period is 4 days, with a fatality rate of 0.00144/day [values based on the reported (official) number of casualties]. The infection fatality rate (IFR) is 0.57 %, and 2.36 % if twice the reported number of casualties is assumed. However, these rates depend on the initially exposed individuals. If approximately nine times more individuals are exposed, there are three times more infected people at the end of the epidemic and IFR = 0.47 %. If we relax these constraints and use a wider range of lower and upper bounds for the incubation and infection periods, we observe that a higher incubation period (13 versus 4.25 days) gives the same IFR (0.6 versus 0.57 %), but nine times more exposed individuals in the first case. Therefore, a precise determination of the fatality rate is subject to the knowledge of the characteristics of the epidemic. We plan to perform again these calculations and publish a short note when the epidemic is over and the complete and precise data is available. Besides the specific example, the analysis proposed in this work shows how isolation measures, social distancing and knowledge of the diffusion conditions help us to understand the dynamics of the epidemic. Hence, the importance to quantify the process to verify the effectiveness of the isolation. There is a long history of mathematical models in epidemiology, going back to the eighteenth century. Bernoulli (1760) used a mathematical method to evaluate the effectiveness of the techniques of variolation against smallpox, with the aim of influencing public health policy. Most of the models are compartmental models, with the population divided into classes and assumptions about the time rate of transfer from one class to another (Hethcote, 2000; Brauer, 2017) . We consider a Susceptible-Exposed-Infected-Removed (SEIR) model to describe the spread of the virus and compute the number of infected and dead individuals. The SEIR model has many versions and the mathematical treatment can be found, for instance, in Hethcote (2000) , Keeling and Rohani (2008) and Diekmann et al. (2013) among others. The goal is to compute the number of infected, recovered and dead individuals on the basis of the number of contacts, probability of the disease transmission, incubation period, recovery rate and fatality rate. The epidemic disease model predicts a peak of infected and dead individuals per day as a function of time, and assumes that births and natural deaths are balanced, since we are dealing with a very short period of time. The population members solely decrease due to the disease dictated by the fatality rate of the disease. The differential equations are solved with a forward-Euler scheme. When no vaccine is available, the isolation of diagnosed infectives and social distancing are the only control measures available. We consider an SEIR epidemic disease model (e.g., Hethcote, 2000; Al-Showaikh and Twizell, 2004; Keeling and Rohani, 2008; Diekmann et al., 2013) . The total (initial) population, N 0 , is categorized in four classes, namely, susceptible, S(t), exposed, E(t), infected, I(t) and recovered, R(t), where t is the time variable. The governing differential equations arė S = Λ − µS − βS I N , where N = S + E + I + R ≤ N 0 in this case, and a dot above a variable denotes time differentiation. Equations (1) are subject to the initial conditions S(0), E(0), I(0) and R(0). The parameters are defined as: All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 April 24, . . https://doi.org/10.1101 April 24, /2020 doi: medRxiv preprint Λ: Per-capita birth rate. µ: Per-capita natural death rate. α: Virus induced average fatality rate. β: Probability of disease transmission per contact (dimensionless) times the number of contacts per unit time. : Rate of progression from exposed to infected (the reciprocal is the incubation period). γ: Recovery rate of infected individuals (the reciprocal is the infection period). having units of (1/T), with T: time. The scheme is illustrated in Figure 1 . The choice Λ = µ = 0 and = ∞ gives the classical SIR model (e.g., d 'Onofrio, 2015) , while if Λ and µ are not zero, the model is termed endemic SIR model (e.g., Allen. 2017 ). However, the SIR model has no latent stage (no exposed individuals) and then it is inappropriate as a model for diseases with an such as that of the COVID-19. Let us clarify better the meaning of each quantity. N is the total number of live humans in the system at time t. S is the number of humans susceptible to be exposed and E is the actual number of exposed individuals (a class in which the disease is latent); people go from S to E depending on the number of contacts with I individuals, multiplied by the probability of infection (β) (see Figure 1 , where βI/N is the average number of contacts with infection per unit time of one susceptible). The other processes taking place at time t are: exposed (E) become infected (I) with a rate and infected recover (R) with a rate γ. Recovered means immune and asymptomatic, but slightly infectious, humans. These individuals do not flow back into the S class, as lifelong immunity is assumed, but it remains to be seen whether the recovered patients from COVID-19 will develop antibodies and achieve lifelong protection. The reciprocals −1 and γ −1 are the average disease incubation and infection periods, respectively. Λ is the rate of birth and µ is the natural rate of death, both per unit time. The reciprocal µ −1 , interpreted as the normal life expectancy (e.g., 70 yr), refers to average normal deaths (e.g., natural deaths, by normal flu, accidents, etc), not related to the infectious disease. These quantities describe a model with vital dynamics (endemic model), which have an inflow of births into the S class at rate Λ and deaths into the other classes at rates µS, µI and µR (see Figure 1 ). If Λ = µN , the deaths balance the newborns and N is constant (but including casualties). The number of live people at time t is N (t) = S(t) + E(t) + I(t) + R(t), that can be lower or higher than N 0 depending on the value of Λ and µ. In this case, it is lower than N 0 . One of the key parameters, besides β, is α that represents the disease-related fatality rate (Chowell et al., 2003; Zhang et al., 2013) . In a very fast pandemic, we may assume that there are no births and normal deaths (or they balance and Λ = µN ), but deaths All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 due to the fatality rate of the disease. This rate is an average, because the model does not take into account the age (a far higher portion of old people die from the disease than young people), the patients preexisting conditions and the healthcare quality. In summary, susceptible persons enter the exposed class with a rate proportional to β and remain there for a mean incubation period −1 , i.e., those already infected with the disease but not able to transmit it are in the exposed class and progress to the infectious class, to recover at the rate γ and die at the rate α. The dead population as a function of time is D(t) = N 0 − N (t), whereas the curve giving the dead people per unit time iṡ (2) Another equivalent approach is an SEIDR model (e.g., De la Sen et al., 2017) , where we have to addḊ to equations (1). In Keeling and Rohani (2008, Section 2.2 where ρ is the per capita probability of dying from the infection. It can easily be shown that equations (2) and (3) are equivalent if births and natural deaths compensate. The basic reproduction ratio, R 0 , is the classical epidemiological measure associated with the reproductive power of the disease. For the SEIR model, it is (Diekmann and Heesterbeek, 2000; Zhang et al., 2013) . It gives the average number of secondary cases of infection generated by an infectious individual. Therefore, it is used to estimate the growth of the virus outbreak. R 0 provides a threshold for the stability of the disease-free equilibrium point. When R 0 < 1, the disease dies out; when R 0 > 1, an epidemic occurs. The behaviour of SEIR models as a function of R 0 can be found, for instance, in Al-Sheikh (2012). The infection fatality rate (IFR) is based on all the population that has been infected, i.e., including the undetected individuals and asymptomatic. In terms of the recovery and fatality rates, we have All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 since the total humans that have been infected is the sum of the recovered and dead individuals, where the subscript refers to the end of the epidemic (t → ∞). It can easily be shown that using the last equation (1) (with µ = 0) and equation (2), we obtain where this relation holds at all times, not only at the end of the epidemic. On the other hand, the case fatality rate (CFR) considers the number of deaths related to the diagnosed individuals, and it is always CFR > IFR, since the number of diagnosed individuals is lower than the denominator of equation (5). The CFR is time dependent and is the usually reported value. We solve the differential equations (1) by using a forward Euler finite-difference scheme (e.g., Carcione, 2014) , discretizing the time variable as t = ndt, where n is a natural number and dt is the time step. Equations (1) and (2) become after discretization: whereḊ n is the number of dead people only in the specific day n. This algorithm yields positive and bounded solutions [e.g., see Brauer (2017) and Problem 1.42(iv) in Diekmann et al. (2013) ], and the system converges to an equilibrium, i.e., S n + R n + D n = S∞ + R∞ + D∞ = N 0 for t → ∞. Let us consider the following base parameters as an example to analyze the results by varying some of them. N 0 = 10 million, α = 0.006/day, β = 0.75/day, γ = (1/8)/day, = (1/3)/day, Λ = µN = 0 (balance of births and natural deaths); and initial conditions: Table 1 ) for SARS and implies an average disease incubation All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 (latent period) of 3 days and an infection period of 8 days. The data correspond to imperfect isolation conditions among individuals and an epidemic situation (high β, R 0 = 5.72 > 1). The time step of the Euler scheme to solve the discretized equations (7) is dt = 0.01 day. Figure 2 shows the number of individuals in the different classes (a), and the total number of dead people (D) and the number of dead people per specific day (Ḋ) (b). As can be seen, the peak of dead individuals per day is reached at day 30. The high values in Figure 2b do not consider home isolation and social distancing measures. Hereafter, we vary the parameters and plot the infected (I) individuals, i.e., excluding those who are incubating the disease (E). Figure 3 shows the number of infected individuals for R 0 > 1 (a) and R 0 ≤ 1 (b), where all the other parameters are kept constant unless β: (µ = 0). We recall here that β is the probability of transmission times the number of contacts per unit time. Basically, reducing β (and R 0 ) the peak decreases in intensity but moves to later times for R 0 higher than 1 (Figure 3a ), although it is wider. There is a significant reduction in the number of infected individuals for R 0 ≤ 1, meaning that home isolation is very effective below a given threshold. The effect of the initially exposed individuals are shown in Figure 4 for two sets of values of R 0 , greater (a) and less (b) than 1. Figure 4a indicates that more exposed people does not mainly affects the intensity of the peak, but anticipates the spread of the epidemic, so that the location of the peak is highly dependent on E(0). On the other hand, Figure 4b shows that for R 0 < 1, the peak location does not change but its intensity does it significantly, with more exposed, more infected. In order the epidemic to start, there should be at least one exposed or one infected individual. Figure 5 indicates that the incubation period (1/ ) has also an impact on the results. If R 0 > 1 (5a), increasing the period from 3 to 9 days decreases the maximum number of infected individuals by almost half and delays the spread of the epidemic, but the peak is wider. If R 0 < 1, the curves behave similarly, but there are much less infected cases. The initially infected individuals (from one to ten thousand) has no apparent effect on the results, as can be seen in Figure 6 , but this is not the case when we deal with the real case history (see next subsection). The effects of the infection period are shown in Figure 7 , where, as expected, increasing this quantity delays the epidemic when R 0 > 1. Below R 0 = 1, the number of infected individuals decreases substantially. Let us assume now that isolation precautions have been imposed and after day 22 R 0 changes from 5.72 to 0.1 [a change of β according to equation (8)], and consider the same parameters to produce Figure 2 . The results are shown in Figure 8 , where the All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10.1101/2020.04.20.20072272 doi: medRxiv preprint peak has moved from day 30 to day 25, with a significant slowing in the number of new cases. The total number of dead individuals has decreased, and the number of dead individuals per day at the peak has decreased from 22 K to 13 K, approximately. Extreme isolation after imperfect isolation anticipates the process. Figure 9 shows the results if the isolation measures start two days before, at day 20 instead of day 22. The number of casualties decreased from 220 K to 155 K. Next, we attempt to model the COVID-19 epidemic in Lombardy (Italy). The time of writing is day 55 (April 19) and the availability of data allows us to perform a relatively reliable fit of the total number casualties from day 1 to date. To predict with high accuracy the behaviour of the epidemic is nearly impossible due to many unknown factors, e.g., the degree of spacial distancing, lack of knowledge of the probability of the disease transmission, characteristics of the disease and parameters of the epidemic. Uncertainties are related to the parameter β that varies with time, while the others are assumed to lie between certain bounds and also contribute to the error. Relative predictions of the trend require an analysis of the data, particularly to define the variation of β and R 0 with time. We do not assume a specific continuous function, but a general approach should consider a partition into discrete periods, , guided by the measures taken by the state and the behaviour of the population. In this case, t 0 = 1 day, t 1 = 22 day and t 2 = 35 day, i.e., L = 3, since after t 1 (March 17), home isolation, social distancing and partial Nation lockdown started to be effective, as indicated by an inflection point in the curve of casualties per day (see below), although it is debatable that the Italian government followed the same rules as in Wuhan, China. We also observe that at t 2 (March 30), the curve starts to bend downwards and reach a "peak". This partition in three periods is valid to date, but the trend can have an unpredictable behavior due to the factors mentioned above, a too early removal of the lockdown conditions, etc. The reported infected people cannot be used for calibration, because these data cannot be trusted. The hospitalization numbers cannot be considered to be representative of the number of infected people and it is largely unknown at present the number of asymptomatic, undiagnosed infections. However, we are aware that even using the number of casualties is uncertain, since there can be an under-ascertainment of deaths, but the figures cannot vary as much as the error related to the infected individuals. Hence, the reported number of deceased people could possibly be underestimated due to undeclared cases. This number depends on the country (quality of the health system) and average age of the population, but it is certain that this novel virus is more deadly All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10. 1101 /2020 and spreads more quickly than seasonal flu. Moreover, authorities make a distinction between a death occurred "with the co-action" of the virus and the death "caused by" the virus. Indeed, only a small percentage of the casualties were in healthy conditions prior the infection and most of the patients were already affected by other illnesses (eg. diabetes, dementia, cancer, stroke). Therefore, we also consider cases where 100 % more people have actually died per day, compared to the official figures. In order to accomplish the fit, we use the simulated annealing algorithm developed by Goffe et al. (1994) . The Fortran code can be found in: https://econwpa.ub. uni-muenchen.de/econ-wp/prog/papers/9406/9406001.txt. The fit is based on the L 2 -norm and yields α, β 1 (before t 0 ), β 2 (after t 0 ), β 3 (after t 1 ), , E(0) and γ from the beginning of the epidemic (day 1, February 25) to date (day 55, April 19), i.e., seven free parameters. We use the total number of deaths for the calibration. Table 1 shows the constraints, initial values and results for different cases, where Cases 1 and 2 correspond to approximately nine times less exposed individuals at the beginning of the epidemic, and Cases 2 and 3 assume double casualties. Cases 4 and 5 consider a wider range of the lower and upper bounds for the incubation and infection periods ( −1 and γ −1 ). The last column is not a variable but indicates the infected individuals at the end of the epidemic, i.e., I∞ = R∞ + D∞ ≈ R∞. Figure 10 shows the curves of Case 1 compared to the data (black dots), with IFR = 0.57 % and R 0 decreasing from 3 to 0.64 at the end of the epidemic. The final infected individuals are 2.37 million people (see Figure 11 and Table 1 ). The inflection point at day 22 (Figure 10b) indicates that the isolation measures started to be effective. Strict isolation could not be achieved at day 22 due to several reasons and there is a reasonable delay of a few days before it can be implemented (day 35). The total number of casualties is approximately 14000 and the effective duration of the epidemic is about 90 days. Recent data reveal that the duration of the Wuhan epidemic was almost 60 days (Wu et al., 2020, Fig . 1b) , a shorter period favoured by the very strict isolation measures applied in that city. Case 2, that considers twice more casualties, and whose results are shown in Figures 12 and 13 , has a high fatality rate, IFR = 2.36 %, and 1.14 million infected people. If the exposed individuals are much higher (Case 3), we obtain IFR = 0.47 % and 5.9 million infected people (see Figures 14 and 15 and Table 1 ). The calculations indicate the uncertainty related to the initially exposed individuals, i.e., those that are incubating the disease. If we modify the constraints and use a wider range of lower and upper bounds, the results are those of Cases 4 and 5 in Table 1 . Case 4 has slightly higher incubation and infection periods compared to Case 1, but a higher IFR (2.25 % versus 0.57 %), whereas more exposed individuals yield an incubation period of 13 days and lower IFR All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10. 1101 /2020 . The results are shown in Figures 16 and 17 (Cases 4 and 5, respectively). Case 6 considers that initially there are a few exposed individuals (we start with one). The algorithm gives a very good fit of the data (plot not shown) with IFR = 3.6 %, comparable periods to Case 1 and 0.63 M infected individuals. Finally, also Case 7, that considers I(0) = 1 and starts with one exposed individual, yields a good fit (plot not shown), but the IFR is too high. These calculations indicate the uncertainty in the determination of the parameters of the epidemic. The values in Table 1 can be compared to figures reported in the literature. The fatality rate and IFR depend on the age of the population. Verity et al. (2020, for China an IFR = 0.657 % but over 60 yr age this rate is 3.28 %. If the number of infected people is several times higher than the reported cases, the fatality rate could be considerably less than the official one, suggesting that this disease is less deadly than SARS and MERS, although much more contagious. Read et al. (2020) report a mean value R 0 = 4, while Wu et al. (2020) obtain values between 1.8 and 2. According to Chowell et al. (2003) , IFR = 4.8 % for SARS, and Verity et al. (2020) state that the average case fatality rate (CFR) of SARS is higher than that of COVID-19, with the latter approximately 1.38 % (their IFR is 0.657 %). However, this virus seems to be much more contagious. The meaning of α −1 is the life expectancy of an individual in the infected class, i.e, if α = 0.00144/day (Case 1), the expectancy is 694 days. There are more complex versions of the SEIR model as, for instance, including a quarantine class and a class of isolated (hospitalized) members (Brauer and Castillo-Chavez, 2012) , or generalizing the diffusion equations (1) with the use of temporal fractional derivatives. The replacement of the first-order temporal derivative by a Caputo fractional derivative of non-natural order provides an additional parameter to fit the data (e.g., Caputo et al., 2011; Mainardi, 2010; Chen et al., 2020) . Furthermore, the model can be made two-dimensional by including the spatial diffusion of the virus (e.g., Naheed et al., 2014) . Moreover, the model can be improved by including others classes. De la Sen et al. (2017) propose an SEIADR model, where A are asymptomatic infectious and D are are dead-infective. In other models, recovered can become again susceptible (e.g., Xia et al., 2016) and, in addition, there are stochastic models (Allen, 2017) , although the calibration becomes extremely difficult with the incomplete data provided by the authorities and the high number of parameters to be found. At the end of the epidemic, more precise information about the parameters will be available and the complete data can be used to evaluate the development of β (and R 0 ) with time. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10.1101/2020.04.20.20072272 doi: medRxiv preprint The outbreak of a pandemic can have catastrophic consequences, not only from the point of view of the casualties, but also economically. Therefore, it is essential to absolutely avoid it by taking the necessary measures at the right time, something that has not been accomplished in Italy and the rest of the world. According to these calculations, the effective measures are social distancing and home isolation, since there is no health system designed for ordinary circumstances that can be prepared for a pandemic, when the infected individuals grows almost exponentially. As can be seen, the pandemic can develop in a few days and the number of casualties can be extremely high if the fatality rate and contagiousness of the disease are high. Only a few days to take action can make a big difference in the prevention of this disaster. The pandemic and its consequences have been predicted in October 2019 by a group of experts (https://www.politico.com/news/magazine/2020/03/07/ coronavirus-epidemic-prediction-policy-advice-121172), but states ignored the fact and transnational nature of the threat, delaying the necessary measures to avoid the disaster, minimizing in many cases the downsides to their own populations and economies. Moreover, in less than three weeks, the virus has overloaded the healthcare system over northern Italy, in particular Lombardy, where the system cannot support this type of emergency and the authorities are not prepared to deal with the epidemic. individual is introduced into a community. It is essential to simulate the process of infection (and death) in advance, to apply adequate control measures and mitigate the risk of virus diffusion. One of the most used mathematical algorithms to describe the diffusion of an epidemic disease is the SEIR model, that we have applied to compute the number of infected, recovered and dead individuals on the basis of the number of contacts, probability of the disease transmission, incubation and infection periods, and disease fatality rate. A first analysis of the results of the model is based on parameters of the SARS disease and we assume that the parameters do not change during the whole epidemic. Reducing the number of contacts, the peak decreases in intensity but moves to later times, although it is wider. Moreover, more exposed people does not affect the intensity of the peak, but anticipates the epidemic. The incubation period has also an impact on the results, with higher values delaying the epidemic. The dependence on the initial number of infected people is weak apparently if R 0 does not change during the epidemic. Increasing the infection period has the opposite effect of increasing the incubation All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10.1101/2020.04.20.20072272 doi: medRxiv preprint period. Moreover, the day when the isolation starts is important, since only two days makes a big difference in the number of casualties. The Lombardy modeling assumes ten million of individuals and has been calibrated on the basis of the total number of casualties. The results show that the peak occurs after 37 days with a final number of dead individuals depending on the reproduction ratio R 0 . With the present available data, this number is almost 14000. Up to day 55 (April 19, the day of writing), the reproduction ratio is 3 before March 17 (day 22), 1.38 [between March 18 and March 30 (day 35)] and 0.64 after March 30, whereas the fatality rate is 0.00144/day (IFR = 0.57 %). We have also doubled the number of casualties and obtained IFR = 2.36 % and 0.47 %, with the second value corresponding to nine times more exposed individuals. These values are obtained by constraining the incubation and infection periods to values reported in the literature. If we relax these constraints and use a wider range of lower and upper bounds, we obtain slightly higher incubation and infection periods compared to the first case, but a much higher IFR (2.25 % versus 0.57 %), while much more exposed individuals yields an incubation period of 13 days and a lower IFR (0.6 %). Therefore, a precise determination of the fatality rate is subject to the knowledge of the parameters of the epidemic and characteristics of the disease, and it is clear from these calculations that the usefulness of simple models to predict is limited, and that their main role is to help in our understanding of the dynamics of the epidemic. We plan to perform again these calculations and publish a short note when the epidemic is over and the complete and precise data is available. Models can be used to predict and understand how an infectious disease spreads in the world and how various factors affect the dynamics. Even if the predictions are inaccurate, it has been clear to scientists from many decades to date that quarantine, social distancing and the adoption of very strict health and safety standards are essential to stop the spreading of the virus. Quarantine was even implemented in medieval times to fight the black death before knowing the existence of viruses. In this sense, this pandemic reveals the failure of policy makers, since it is well known from basic modeling results that anticipating those measures can save thousand of lives and even prevent the pandemic. The interface of science, society and politics is still uneasy, even in highly developed countries, revealing a disregard for scientific evidence. Moreover, one of the consequences is that some of these countries do not invest sufficiently in R&D and must acquire the new technology overseas at a much higher cost. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 verified the figures and contributed to the analysis; CB aided in the solution of the differential equations; JB contributed to the discussion. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Acknowledgements: JMC dedicates this work to Antonio Buonomo, who passed away of COVID-19 on March 21, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 understanding infectious disease dynamics. Princeton Series in Theoretical and Computational Biology. Princeton University Press, Princeton. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 Manuscript contribution to the field: We implement an SEIR model to evaluate the date of the infection peak of the COVID-19 epidemic, that includes the disease fatality rate to estimate the number of casualties per day. To our knowledge, this is the first time that this model is calibrated with the number of casualties. The simulation attempts to provide a simple procedure to model the coronavirus diffusion in a given region. The example regards the epidemic in the Lombardy province (Italy) which is taking place at the time of this writing, but can be applied in general. We show how the date of the peak and the number of casualties are affected by the effectiveness of the home isolation, incubation and infection periods, probability of transmission and initially exposed individuals. The results from the procedure, that intends to be the basis for a real case study, where it is clear the effect of each parameter and variable in the dynamic evolution of the epidemic. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. . https://doi.org/10.1101/2020.04.20.20072272 doi: medRxiv preprint Table 1 . Constraints and initial-final values of the inversion algorithm. Variable Lower bound 10 −5 0.5 10 −6 10 −6 3 3 10 3 Upper bound 10 −1 0.9 10 3 10 3 6 6 2 × 10 5 Initial value 0.006 0.75 0.5 0.5 5 5 10 4 /10 5 1. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Fig. 1 A typical SEIR model. The total population, N , is categorized in four classes, namely, susceptible, S, exposed E, infected I and recovered R (e.g., Chitnis et al., 2008) . Λ and µ correspond to births and natural deaths independent of the disease, and α is the fatality rate. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 Fig. 4 Infected individuals for different values of the initially exposed individuals, corresponding to R 0 greater (a) and smaller (b) than 1. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Table 1 . The peak can be observed at day 37 (April 1). All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Table 1 . All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Table 1 . All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Table 1 ). All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020. Table 1 ). All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted April 24, 2020 . . https://doi.org/10.1101 /2020 Modeling and analysis of an SEIR epidemic nodel with a limited resource for treatment Grand challenge in human/animal virology: Unseen, smallest replicative entities shape the whole globe A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis One-dimensional measles dynamics Essai d'une nouvelle analyse de la mortalité causée par la petite vérole et des avantages de l'inoculation pour la prévenir, Mémoires de Mathématiques et de Physique Mathematical epidemiology: Past, present, and future Mathematical Models in Population Biology and Epidemiology Wave simulation in biological media based on the Kelvin-Voigt fractional-derivative stress-strain relation Wave Fields in Real Media. Theory and numerical simulation of wave propagation in anisotropic, anelastic, porous and electromagnetic media The reconstruction and prediction algorithm of the fractional TDD for the local outbreak of COVID-19 Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model SARS outbreak in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism Mathematical tools for On a new epidemic model with asymptomatic and dead-infective subpopulations with feedback controls useful for Ebola disease Dynamic behaviour of a discrete-time SIR model with information dependent vaccine uptake Global optimization of statistical functions with simulated annealing The mathematics of infectious diseases Modeling infectious diseases in humans and animals The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application Fractional Calculus and Waves in Linear Viscoelasticity Numerical study of SARS epidemic model with the inclusion of diffusion in the system Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions Pale rider: The spanish flu of 1918 and how it changed the world Estimates of the severity of coronavirus disease 2019: a model-based analysis Dynamics of a delayed SEIQ epidemic Model Estimating clinical