Global analysis of the COVID-19 pandemic using simple epidemiological models
aa r X i v : . [ q - b i o . P E ] M a y Global analysis of the COVID-19 pandemic using simple epidemiological models
Jos´e Enrique Amaro, ∗ J´er´emie Dudouet, † and Jos´e Nicol´as Orce ‡ Departamento de F´ısica At´omica, Molecular y Nuclear and Instituto Carlos I de F´ısica Te´orica y Computacional,Universidad de Granada, E-18071 Granada, Spain. Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3,IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France Department of Physics & Astronomy, University of the Western Cape, P/B X17 Bellville ZA-7535, South Africa. (Dated: May 15, 2020)Several analytical models have been used in this work to describe the evolution of death casesarising from coronavirus (COVID-19). The Death or ‘D’ model is a simplified version of the SIR(susceptible-infected-recovered) model, which assumes no recovery over time, and allows for thetransmission-dynamics equations to be solved analytically. The D-model can be extended to describevarious focuses of infection, which may account for the original pandemic (D1), the lockdown (D2)and other effects (Dn). The evolution of the COVID-19 pandemic in several countries (China,Spain, Italy, France, UK, Iran, USA and Germany) shows a similar behavior in concord with theD-model trend, characterized by a rapid increase of death cases followed by a slow decline, which areaffected by the earliness and efficiency of the lockdown effect. These results are in agreement withmore accurate calculations using the extended SIR model with a parametrized solution and moresophisticated Monte Carlo grid simulations, which predict similar trends and indicate a commonevolution of the pandemic with universal parameters.
PACS numbers:Keywords: COVID-19, death model,ESIR model,Monte Carlo Planck model
I. MOTIVATION
The SIR (susceptible-infected-recovered) model iswidely used as first-order approximation to viral spread-ing of contagious epidemics [1], mass immunization plan-ning [2, 3], marketing, informatics and social networks [4].Its cornerstone is the so-called “mass-action” principleintroduced by Hamer, which assumes that the course ofan epidemic depends on the rate of contact between sus-ceptible and infected individuals [5]. This idea was ex-tended to a continuous time framework by Ross in hispioneering work on malaria transmission dynamics [6–8], and finally put into its classic mathematical form byKermack and McKendric [9]. The SIR model was furtherdeveloped by Kendall, who provided a spatial generaliza-tion of the Kermack and McKendrick model in a closedpopulation [10] (i.e. neglecting the effects of spatial mi-gration), and Bartlett, who – after investigating the con-nection between the periodicity of measles epidemics andcommunity size – predicted a traveling wave of infectionmoving out from the initial source of infection [11, 12].More recent implementations have considered the typicalincubation period of the disease and the spatial migrationof the population.The
COVID-19 pandemic has ignited the submission ofmultiple manuscripts in the last weeks. Most statisticaldistributions used to estimate disease occurrence are ofthe binomial, Poisson, Gaussian, Fermi or exponential ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] types. Despite their intrinsic differences, these distribu-tions generally lead to similar results, assuming indepen-dence and homogeneity of disease risks [13].In this work, we propose a simple and easy-to-use epi-demiological model – the Death or D model [14] – thatcan be compared with data in order to investigate theevolution of the infection and deviations from the pre-dicted trends. The D model is a simplified version ofthe SIR model with analytical solutions under the as-sumption of no recovery – at least during the time of thepandemic. We apply it globally to countries where theinfestation of the
COVID-19 coronavirus has widespreadand caused thousands of deaths [15, 16].Additionally, D-model calculations are benchmarkedwith more sophisticated and reliable calculations usingthe extended SIR (ESIR) and Monte Carlo Planck(MCP) models – also developed in this work – whichprovide similar results, but allow for a more coherentspatial-time disentanglement of the various effectspresent during a pandemic. A similar ESIR model hasrecently been proposed by Squillante and collaboratorsfor infected individuals as a function of time, based onthe Ising model – which describes ferromagnetism in sta-tistical mechanics – and a Fermi-Dirac distribution [17].This model also reproduces a posteriori the
COVID-19 data for infestations in China as well as other pandemicssuch as Ebola, SARS, and influenza A/H1N1.The SIR model considers the three possible states ofthe members of a closed population affected by a conta-gious disease. It is, therefore, characterized by a systemof three coupled non-linear ordinary differential equa-tions [18], which involve three time-dependent functions: • Susceptible individuals, S ( t ), at risk of becominginfected by the disease. • Infected individuals, I ( t ). • Recovered or removed individuals, R ( t ), who wereinfected and may have developed an immunity sys-tem or die.The SIR model describes well a viral disease, whereindividuals typically go from the susceptible class S tothe infected class I , and finally to the removed class R .Recovered individuals cannot go back to be susceptible orinfected classes, as it is, potentially, the case of bacterialinfection. The resulting transmission-dynamics systemfor a closed population is described by dSdt = − λSI, (1) dIdt = λSI − βI, (2) dRdt = βI, (3) N = S ( t ) + I ( t ) + R ( t ) , (4)where λ > β > N is the fixed population size,which implies that the model neglects the effects of spa-tial migration. Currently, there is no vaccination avail-able for COVID-19 , and the only way to reduce the trans-mission or infection rate λ – which is often referred to as“flattening the curve”– is by implementing strong socialdistancing and hygiene measures.The system is reduced to a first-order differential equa-tion, which does not possess an explicit solution, butcan be solved numerically. The SIR model can then beparametrized using actual infection data to solve I ( t ), inorder to investigate the evolution of the disease. In theD model, we make the drastic assumption of no recov-ery in order to obtain an analytical formula to describe –instead of infestations – the death evolution by COVID-19 . This can be useful as a fast method to foresee theglobal behavior as a first approach, before applying moresophisticated methods. We shall see that the resultingD model describes well enough the data of the currentpandemics in different countries.
II. THE DEATH OR D MODEL
The main assumption of the D model is the absence ofrecovery from coronavirus, i.e. R ( t ) = 0, at least duringthe pandemic time interval. This assumption may bereasonable if the spreading time of the pandemic is muchfaster than the recovery time, i.e. λ ≫ β . The SIRequations are then reduced to the single equation of thewell-known SI model, dIdt = λ ( N − I ( t )) I ( t ) , (5) which represents the simplest mathematical form of alldisease models, where the infection rate is proportionalto both the infected, I , and susceptible individuals N − I .Equation 5 is trivially solved by multiplying by dt anddividing by ( N − I ) I , dI ( N − I ) I = λdt, (6)or (cid:18) N − I + 1 I (cid:19) dI = λN dt. (7)Integrating over an initial t = 0 and final t we obtain ℓn I ( t ) N − I ( t ) − ℓn I N − I = λN ( t − t ) , (8)where I = I ( t ). Taking the exponential on both sides I ( t ) N − I ( t ) = I N − I e λN ( t − t ) . (9)Finally, solving this algebraic equation we obtain the so-lution I ( t ) I ( t ) = N I e λN ( t − t ) N − I + I e λN ( t − t ) , (10)which can be written in the form I ( t ) = I e ( t − t ) /b − C + C e ( t − t ) /b , (11)where we have defined the constants b = 1 λN , C = I N . (12)The parameter b is the characteristic evolution time ofthe initial exponential increase of the pandemic. Theconstant C is the initial infestation rate with respect tothe total population N . Assuming C ≪
1, Eq. 11 yields I ( t ) = I e ( t − t ) /b C e ( t − t ) /b . (13)In order to predict the number of deaths in the D modelwe assume that the number of deaths at some time t isproportional to the infestation at some former time τ ,that is, D ( t ) = µI ( t − τ ) , (14)where µ is the death rate, and τ is the death time. Withthis assumption we can finally write the D-model equa-tion as D ( t ) = a e ( t − t ) /b c e ( t − t ) /b , (15)where a = µI e − τ/b , c = C e − τ/b , and a/c yields thetotal number of deaths predicted by the model. This is a = 70 . b = 6 .
35 days, c = 0 . D E A T H S CHINA D ( t ) = a exp( t/b )1 + c exp( t/b ) 60504030201003500300025002000150010005000
10 20 30 40 50 60 70
Days D ea t h s / da y datafit January 22 a exp(t/b)b(1+c exp(t/b)) February 15
D’(t) = a = 76.9666 b = 6.4094 days c = 0.0235
FIG. 1: Fits to total (left panel) and daily (right panel) deaths by COVID-19 in China using the D ( t ) and D ′ (t), respectively.The dashed curve shows a fit to the daily deaths using the parameters determined to fit the total deaths (top of left panel),which provides similar results to an independent fit (parameters on the top right), given the statistical fluctuations in the dailyrates. Data are taken from [19]. the final equation for the D-model, which presents a sim-ilar shape to the well-known Woods-Saxon potential forthe nucleons inside the atomic nucleus or the bacterialgrowth curve. The rest of the parameters, µ , τ , I and N are embedded in the parameters a, b, c , which repre-sent space-time averages and can be fitted to the timelyavailable data.In Fig. 1, we present the fit of the D-model to the COVID-19 death data for China, where its evolutionhas apparently been controlled and the D function hasreached the plateau zone, with few increments over time,or fluctuations that are beyond the model assumptions.This plot shows the duration of the pandemic – abouttwo months to reach the top end of the curve – and theagreement, despite the crude assumptions, between dataand the evolution trend described by the D-model. Thisagreement encourages the application of the D modelto other countries in order to investigate the differenttrends.
A. Evolution of D-model parameters
In order to get insight into the stability and uncertaintyof our predictions, Fig. II shows the evolution of a , b , and c and other model predictions from fits to the daily datain Spain. The meaning of these quantities is explainedbelow: • The parameter a is the theoretical number ofdeaths at the day corresponding to t = 0. In gen-eral, it differs from the experimental value and canbe interpreted as the expected value of deaths thatday. Note that experimental data may be subjectto unknown systematic errors and different count-ing methods. • The parameter b , as mentioned above, is the char- DAYS a Deaths on day 1 3028262422201816504540353025 DAYS b [ d a y s ] Evolution time 30282624222018164 . . . . . . c Inverse death factor 30282624222018160 . . . . . . . .
002 DAYS a / c Total deaths 3028262422201816160001500014000130001200011000DAYS T [ d a y s ] Time to 99% 302826242220181642414039383736 DAYS T [ d a y s ] Time to 95% 302826242220181636353433323130
FIG. 2: Evolution of a, b and c parameters and various pre-dictions of the D-model as a function of time (days). acteristic evolution time. During the initial expo-nential behavior, it indicates the number of daysfor the number of deaths to double. Moreover, 1 /b is proportional to the slope of the almost linear be-havior in the mid region of the D function. Thatbehavior can be obtained by doing a Taylor expan-sion around t = − b ℓn c and is given by D ( t ) ≃ c (cid:18) − ℓn c (cid:19) + t bc . (16) • The parameter c is called the inverse dead factor because D ( t → ∞ ) = a/c provides the asymptoticor expected total number of deaths. • The times T and T correspond to D =0 . D ( ∞ ) and D = 0 . D ( ∞ ), respectively. Thesetimes are obtained by solving the equation D ( t ) = γa/c , where γ = 0 .
95 or 0.99. The solution of thatequation is t = b ℓn (cid:18) c γ − γ (cid:19) . (17)Figure II shows the stable trend of the parameters be-tween days 19 to 24 (corresponding to March 27–30),right before reaching the peak of deaths cases, which oc-curred in Spain around April 1. Such stability validatesthe D-model predictions during this time. However, arapid change of the parameters is observed, especiallyfor a , once the peak is reached, drastically changing theprediction of the number of deaths given by a/c . Thissudden change results in the slowing down of deaths perday and longer time predictions T and T .The parameters of the D model correspond to averagevalues over time of the interaction coefficients betweenindividuals, i.e. they are sensitive to an additional exter-nal effect on the pandemic evolution. These may includethe lockdown effect imposed in Spain in March 14 andother effects such as new sources of infection or a sud-den increase of the total susceptible individuals due tosocial migration and large mass gatherings [20]. It is notpossible to identify a specific cause because its effectsare blurred by the stochastic evolution of the pandemic,which is why any reliable forecast presents large errors. B. The D ′ model One can also determine deaths/day rates by applyingthe first derivative to Eq. 15, D ′ ( t ) = a e ( t − t ) /b b (1 + c e ( t − t ) /b ) , (18)which allows for a determination of the pandemics peakand evolution after its turning point. The D model de-scribes well the cumulative deaths because the sum ofdiscrete data reduce the fluctuations, in the same way as the integral of a discontinuous function is a continu-ous function. However, the daily data required for D ′ have large fluctuations – both statistical and systematic– which normally gives a slightly different set of param-eters when compared with the D model. a = 84 . b = 4 .
61 days, c = 0 . D E A T H S / D AY
16 April4 AprilSPAIN a exp( t/b ) b (1 + c exp( t/b )) a = 93 . b = 4 .
61 days, c = 0 . D E A T H S / D AY
16 April4 AprilSPAIN D ( t ) − D ( t −
1) 50403020100120010008006004002000
FIG. 3: Predictions of the D model for D ′ ( t ) in Spain accord-ing with the data collected up to April 5. Using the D model fitted to cumulative deaths allowsto compute deaths/day as D ( t ) − D ( t − ∆ t ) ≃ D ′ ( t )∆ t, (19)where ∆ t = 1 day. Figure 3 shows that Eqs. 18 and 19yield similar parameters, as the time increment is smallenough compared with the time evolution of the D ( t )function. Hence, the first derivative D ′ ( t ) can be usedto describe deaths per day. In addition, Fig. 4 showsthat the parameters may be different for both D and D ′ functions using cumulative and daily deaths, respectively,as shown for Spain on April 5. It is also important tonote that b is directly proportional to the full width athalf maximum ( F W HM ) of the D ′ ( t ) distribution, F W HM = 2 b ℓn (3 + 2 √ ≈ . b. (20)As shown below, the b parameter presents typical val-ues between 4 and 10 for most countries undergoing theinitial exponential phase, which yields a minimum andmaximum time of 14 and 35 days, respectively, betweenthe two extreme values of the F W HM . C. D n model with two or more channels of infection Some models [21] include changes in the transmissionrate due to various interventions implemented to containthe outbreak. The simple D model does not allow to dothis explicitly, but changes in the spread can be takeninto account by considering the total D or D n functionas the sum of two or more independent D-functions withdifferent parameters, which may reveal the existence ofseveral independent sources, or virus channels. An exam-ple is shown in Fig. 5, where the two-channel function D ′ = D ′ ( a, b, c ) + D ′ ( a , b , c ) , (21) a = 49 . b = 4 . c = 0 . D E A T H S
16 AprilSPAIN, 5 April a exp( t/b )1 + c exp( t/b ) 40353025201510501600014000120001000080006000400020000 a = 78 . b = 4 .
53 days, c = 0 . D E A T H S / D AY
16 April5 AprilSPAIN a exp( t/b ) b (1 + c exp( t/b )) FIG. 4: Comparison of parameters fitted to D ( t ) and D ′ ( t )in Spain according with the data on April 5. has been fitted with six parameters to the Spanish dataup to April 13. The fit reveals a second, smaller deathpeak, which substantially increase the number of deathsper day and the duration of the pandemic. This is equiv-alent to add a second, independent, source of infectionseveral weeks after the initial pandemic. The second peakmay as well represent a second pandemic phase drivingthe effects of quarantine during the descendant part ofthe curve. a = 35, b = 3 . c = 0 . a = 14, b = 5, c = 0 . DAYS D E A T H S / D AY
20 April13 AprilSPAIN D ′ Model 605040302010010008006004002000
FIG. 5: Predictions of the D ′ model in Spain using a sum oftwo D ′ -functions for data collected up to April 13. Additionally, the cumulative D-function can also becomputed with a two-channel function, D = D ( a, b, c ) + D ( a , b , c ) , (22)which provides, as shown in Fig. 6, a more accurateprediction for the total number of deaths and clearly il-lustrates the separate effect of both source peaks. It isinteresting to note that for large t , a ≈ a , c ≈ c and b ≈ b . In such a case, the total number of deaths ex-pected during the pandemic is given by D ( ∞ ) = 2 a/c . a = 26, b = 3 . c = 0 . a = 9 . b = 4 . c = 0 . DAYS T O T A L D E A T H S March 30 April 13April 20SPAIN D model 605040302010020000150001000050000 FIG. 6: Predictions of the D model in Spain using a sum oftwo D-functions for data collected up to April 13. D. Estimation of the infected function I ( t ) The D-model can also be used to estimate I ( t ) usingthe initial values of I = I (0) and the total number ofsusceptible people N = S (0). The initial value of N isunknown, and not necessarily equal to the population ofthe whole country since the pandemic started in local-ized areas. Here, we shall assume N = 10 , althoughplausible values of N can be tens of millions. Note thatthe no-recovery assumption of the D model is unrealistic,and this calculation only provides an estimation of thenumber of individuals that were infected at some time,independently of whether they recovered or not.From the definition of D ( t ) in Eq. 14, the followingrelations between the several parameters of the modelwere extracted a = µI e − τ/b , (23) c = I N e − τ/b , (24) b = 1 λN . (25)Solving the first two equations for µ and I we obtain I = N c e τ/b , (26) µ = aN c . (27)Hence, µ can be computed by knowing N . However, toobtain I one needs to know the death time τ . This hasbeen estimated to be about 15 to 20 days for COVID-19 cases, which can be used to compute two estimates of I ( t ). These are given in Fig. 7 for the case of Spain.Since there is no recovery in the D model, the totalnumber of infected people is I ∼ N for large t , i.e. N =10 in our case. In Fig. 7, we have labeled the beginningof the lockdown in Spain (March 15). For τ = 15 days,most of the susceptible individuals were already infectedon that date, and even more for τ = 20 days, as the τ = 20 d τ = 15 d a = 49 . b = 4 . c = 0 . C A S E S InfectedDeaths
March 15 April 6February 1
SPAIN, N = 10 − − − − × × × FIG. 7: Predictions of the D model for the infected function I ( t ) in Spain according to data collected up to April 6. τ = 20 d τ = 15 d a = 49 . b = 4 . c = 0 . D ( t ) / I ( t ) SPAIN, N = 10 . . . . . . . . τ = 20 d τ = 15 d a = 49 . b = 4 . c = 0 . I ( t ) / N April 6March 7
SPAIN 403020100 − − . . . . FIG. 8: (Top panel) Results for D ( t ) /I ( t ) (deaths over in-fected) and (bottom panel) I ( t ) /N (infected over susceptible)according to the D-model. Results are for data collected inSpain up to April 6, assuming N = 10 . pandemic had started almost two months earlier. Mostof the individuals got infected, even if a great part ofthem – approximately 99% – had no symptoms of illnessor disease.Moreover, the top panel of Fig. 8 shows the ratio D ( t ) /I ( t ) (deaths over infected), as given by Eqs. 13 and 15, D ( t ) I ( t ) = aN c e τ/b c e ( t + τ ) /b c e τ/b , (28)which also depends on N and τ . For N = 10 , the ratio D/I increases similarly to the separate functions D and I between the initial and final values, D (0) I (0) = aN c e τ/b , (29) D ( ∞ ) I ( ∞ ) = aN c . (30)These results depend on the total susceptible popula-tion N . However, the ratio of infected with respect tosusceptibles, I/N , is independent on N . This functiondepends only on τ and is shown in the bottom panel ofFig. 8 for τ = 15 and 20 days, which reveals the rapidspread of the pandemic. Accordingly, between 10% and30% of the susceptibles were infected in March 7, and onemonth later (April 6), when the fit was made, all suscepti-bles had been infected. This does not means that the fullpopulation of the country got infected, since the number N is unknown and, for instance, excludes individuals inisolated regions, and it may additionally change becauseof spatial migration, not considered in the model. III. THE EXTENDED SIR MODEL
D-model predictions can be compared with more re-alistic results given by the complete SIR model [9, 11],which is characterized by Eqs. 1, 2, 3 and 4 with initialconditions R (0) = 0, I (0) = I , S (0) = N − I . TheSIR system of dynamical equations can be reduced to anon-linear differential equation. First, dividing Eq. 1 byEq. 3 one obtains, dSdR = − λβ S, (31)which yields the following exponential relation betweenthe susceptible and the removed functions, S = S e − λR/β . (32)Moreover, Eq. 4 provides a relation between the infectedand the removed functions, I = N − S − R = N − S e − λR/β − R, (33)which yields, by inserting into Eq. 3, the final SIR differ-ential equation dRdt = β (cid:16) N − S e − λR/β − R (cid:17) . (34)In order to obtain R ( t ) we only need to solve thisfirst-order differential equation with the initial condition Extended SIR model - no end conditionDAYS D E A T H S / D AY March 30 April 14April 20SPAIN a = 0 . b = 4 . c = 0 . a = 982 b = 0 . c = 1 Extended SIR model D ′ (100) = 10DAYS D E A T H S / D AY March 30 April 15April 20SPAIN a = 0 . b = 4 . c = 0 . a = 881 b = 0 . c = 1 . D ′ (100) = 5DAYS D E A T H S / D AY March 30 April 15April 20SPAIN a = 0 . b = 4 . c = 0 . a = 881 b = 0 . c = 1 . FIG. 9: Fit of the ESIR model to daily deaths in Spain up to April 15 using no boundary condition for the final number ofdeaths (left panel), and with boundary conditions of D ′ (100) = 10 (middle panel) and D ′ (100) = 5 (right panel) deaths/day. R (0) = 0. Moreover, if we normalize the functions S , I and R to 1, S = sN, (35) I = iN, (36) R = rN, (37)so that s + i + r = 1, then r ( t ) verifies drdt = β (cid:16) − s e − λNr/β − r (cid:17) , (38)which can be solved numerically, or by approximatemethods in some cases. In Ref. [9], a solution was foundfor small values of the exponent λN r/β . For the coro-navirus pandemic, however, this number is expected toincrease and be close to one at the pandemic end.At this point, we propose a modification of the stan-dard SIR model. Instead of solving Eq. 38 numericallyand fitting the parameters to data, the solution can beparametrized as r ( t ) = ac + e − t/b , (39)which presents the same functional form as the D-modeland, conveniently, provides a faster way to fit the modelparameters by avoiding the numerical problem of solvingEq. 38. In fact, numerical solutions of the SIR modelpresent a similar step function for R ( t ). Additionally,one can assume that D ( t ) is proportional to R ( t ), andcan also be written as dDdt = a (cid:16) − c e − r/b − r ( t ) (cid:17) , (40)where a , c = s and b = β/ ( λN ) are unknown param-eters to be fitted to deaths-per-day data, together withthe three parameters of the r ( t )-function: a , b , c .Figure 9 shows fits of the ESIR model to daily deathsin Spain during the coronavirus spread. The use of no boundary condition for the number of deaths (left panel)is not an exact solution of the SIR differential equation.A way to solve this problem is to impose the condition D ′ ( ∞ ) = 0, as the number of deaths must stop at sometime. Numerically, it is enough to choose a small valueof D ′ ( t ) for an arbitrary large t . The middle and rightpanels of Fig. 9 show different boundary conditions of D ′ (100) = 10 and D ′ (100) = 5, respectively, which yieldthe same results and the expected behavior for a viraldisease spreading and declining.It is also consistently observed (e.g. see middle andright panels of Fig. 9), that at large t , r ( t ) → ac ≈ N recovers, as we previously inferred from the D model. This, together with the fact that c can alwaysbe adjusted to 1, leaves the ESIR model with essentially4 free parameters to fit to the daily death data; i.e. thesame number of parameters than the original SIR model.As shown in Fig. 10, ESIR fits reproduce well the longflattening behavior observed in UK, USA, Germany orIran, whereas it fails to reproduce the more-pronounceddouble-peak structure typically observed in countries likeFrance, Italy, Spain or Belgium.As previously done with the D model, one can alsoexpand the ESIR model to accommodate this apparentfailure to take lockdown effects into account. Similarly,the ESIR2 model is proposed as,ESIR2( t ) = a (cid:16) − c e − r/b − r ( t ) (cid:17) , (41)with r ( t ) = ac + e − t/b + a ′ c ′ + e − t/b ′ = a a + e − t/b + a a + e − t/b ′ , (42)where we have assumed that a = a ′ and c = 2 a to ac-commodate that r ( ∞ ) → c = 1. Hence, we areleft with five free parameters. FIG. 10: Fit to the data (average of 7 consecutive days up to May 8) of the ESIR and D ′ models in the United Kingdom (left)and France (right). Finally, Fig. 11 shows the comparison between theESIR2 and D ′ fits to real data for countries where COVID-19 has widely spread: Belgium, USA, France,Germany, Iran, Italy, Spain and UK, USA. Death dataare taken from Refs. [19, 22, 23] and consider 7-day aver-age smoothing to correct for anomalies in data collectionsuch as the typical weekend staggering observed in var-ious countries, where weekend data are counted at thebeginning of the next week. Real error intervals are ex-tracted from the correlation matrix. As discussed in Sec-tion 2.3, the reduced D ′ model has been used with a = a and c = c . Although arising from different assumptions,both models provide similar data descriptions and pre-dictions, with slightly better values of χ per degree offreedom for the ESIR2 model. It is also interesting tonote that the reduced ESIR2 model with five parametersyields similar results to the full ESIR2 model, with eightparameters.As data become available, daily predictions vary forboth ESIR2 and the D ′ models. This is because themodel parameters are actually statistical averages overspace-time of the properties of the complex system. Nomodel is able to predict changes over time of these prop-erties if the physical causes of these changes are not in-cluded. The values of the model parameters are only welldefined when the disease spread is coming to an end andtime changes in the parameters have little influence. IV. DISCUSSION OF GLOBAL RESULTS
More sophisticated calculations can be compared withESIR2 and D ′ predictions. In particular, Monte Carlo(MC) simulations have also been performed in this workfor the Spanish case [24], which consist of a lattice of cellsthat can be in four different states: susceptible, infected,recovered or death. An infected cell can transmit thedisease to any other susceptible cell within some randomrange R . The transmission mechanism follows principlesof nuclear physics for the interaction of a particle with a target. Each infected particle interacts a number n of times over the interaction region, according to its en-ergy. The number of interactions is proportional to theinteraction cross section σ and to the target surface den-sity ρ . The discrete energy follows a Planck distributionlaw depending on the ’temperature’ of the system. Forany interaction, an infection probability is applied. Fi-nally, time-dependent recovery and death probabilitiesare also applied. The resulting virus spread for differentsets of parameters can be adjusted from COVID-19 pan-demic data. In addition, parameters can be made timedependent in order to investigate, for instance, the effectof an early lockdown or large mass gatherings at the riseof the pandemic.As shown in Fig. 12, our MC simulations present sim-ilar results to the D ′ model, which validates the use ofthe simple D-model as a first-order approximation. Moredetails on the MC simulation will be presented in a sep-arate manuscript [24]. Interestingly, MC simulations fol-low the data trend up to May 11 without any changes inthe parameters for nearly two weeks. An app for Androiddevices, where the Monte Carlo Planck model has beenimplemented to visualize the simulation is available fromRef. [25].In order to investigate the universality of the pan-demic, it is interesting to compare all countries by plot-ting the D model in terms of the variable ( t − t ) /b ,where t is the maximum of the daily curve given by t max = − b ℓn ( c ). By shifting Eq. 15 by t max = − b ℓn ( c )and dividing by t max = a/c , the normalized D functionis given by, D norm ( t ) = c e ( t − t max ) /b c e ( t − t max ) /b . (43)The left of Fig. 13 shows similar trends for the nor-malized D curves of different countries, which suggests auniversal behavior of the COVID-19 pandemic. Only Iranseems to slightly deviate from the global trend, whichmay indicate an early and more effective initial lockdown.A similar approach can be done for the daily data using
FIG. 11: Daily deaths fitted with the ESIR2 and D ′ models with a cut off of May 8 2020. the D ′ and ESIR2 models, as shown in the middle andright panels of Fig. 13, respectively. Although differentcountries show similar trends, statistical fluctuations inthe daily data do not result in a nice universal behavior as compared with D norm . However, the D ′ and ESIR2plots show that an effective lockdown is characterized byflatter and broader peaks, best characterized the Iraniancase, whereas Spain and Germany present the sharper0 SpainMonte CarloD2 model a = 44, b = 4 . c = 0 . a = 58, b = 7 . c = 0 . DAYS T O T A L D E A T H S March 30 May 11 100806040200300002500020000150001000050000
Monte Carlo P-Model with lockdownDAYS D E A T H S / D AY Monte Carlo19 April
SPAIN 605040302010010008006004002000
FIG. 12: Predictions of the MC and D models in Spain up to May 11. peaks.
V. FINAL REMARKS
The global models considered in this work present somedifferences with respect to other existing models. First,in this work we have tried to keep the models as simpleas possible. This allows to use theoretical-inspired ana-lytical expressions or semi-empirical formulae to performthe data analysis. The use of semi-empirical expressionsfor describing physical phenomena is recurrent in physics.One of the most famous is the semi-empirical mass for-mula from nuclear physics. Of course the free parametersneed to be fitted from known data, but this allowed toobtain predictions for unknown elements.In our case we were inspired by the well known sta-tistical SIR-kind models slightly modified to obtain an-alytical expressions that carry the leading time depen-dence. We have found that the D and D models allowa fast and efficient analysis of the pandemics in the ini-tial and advanced stages. Our results show that the timedependence of the pandemic parameters due to the lock-down can be effectively simulated by the sum of two D-functions with different widths and heights and centeredat different times. The distance between the maximaof the two D-functions should be a measure of the timebetween the effective pandemic beginning and lockdown.In the Spanish case this is about 20 days. Taking intoaccount that lockdown started in March 14, this marksthe pandemic starting time as about February 22. Hadthe lockdown started on that date, the deaths would hadbeen highly reduced. The smooth blending between thetwo peaks provides a transition between the two statisti-cal regimes (or physical phases) with and without lock-down.The Monte Carlo simulation results are in agreementwith our previous analysis with the D and D models.The Monte Carlo generates events in a population of in- dividuals in a lattice or grid of cells. We simulate themovement of individuals outside of the cells and inter-actions with the susceptible individuals within a finiterange. The randon events follow statistical distributionsbased on the exponential laws of statistical mechanicsfor a system of interacting particles, driven by macro-scopic magnitudes as the temperature, and interactionprobabilities between individuals, that can be related tointeraction cross sections.The Monte Carlo simulation spread the virus in space-time, and allows also space-time dependence on the pa-rameters. In this work we have made the simplest as-sumptions, only allowing for a lockdown effect by reduc-ing the range of the interaction starting on a fixed day.This simple modification allowed to reproduce nicely theSpanish death-per-day curve. The lockdown producesa relatively long broadening of the curve and a slow de-cay. Similar MC calculations can be performed in severalcountries to infer the devastating effect of a late lockdownas compared with early lockdown measures. The later isthe case of South Africa and other countries, which havenot reached the exponential growth.The Death and extended SIR models are simple enoughto provide fast estimations of pandemic evolution by fit-ting spatial-time average parameters, and present a goodfirst-order approximation to understand secondary effectsduring the pandemic, such as lockdown and populationmigrations, which may help to control the disease. Sim-ilar models are available [17, 26], but challenges in epi-demiological modeling remain [27–30]. This is a verycomplex system, which involves many degrees of freedomand millions of people, and even assuming consistent dis-ease reporting - which is rarely the case – there remainsan important open question: Can any model predict theevolution of an epidemic from partial data? Or similarly,Is it possible, at any given time and data, to measure thevalidity of an epidemic growth curve? We finally hopethat we have added new insightful ideas with the Death,the extended SIR and Monte Carlo models, which can1 FIG. 13: Universality of the normalized D (left), D ′ (middle) and ESIR2 (right) models. now be applied to any country which has followed theinitial exponential pandemic growth. Acknowledgments
The authors thank useful comments from EmmanuelCl´ement, Araceli Lopez-Martens, David Jenkins, Ra- mon Wyss, Liam Gaffney and Hans Fynbo. Thiswork was supported by the Spanish Ministerio deEconom´ıa y Competitividad and European FEDERfunds (grant FIS2017-85053-C2-1-P), Junta de Andaluc´ıa(grant FQM-225) and the South African National Re-search Foundation (NRF) under Grant 93500. [1] R. M. Anderson, Discussion: the Kermack-McKendrickepidemic threshold theorem. Bulletin of mathematical bi-ology, : 132 (1991).[2] S. Chauhan1, O. P. Misra and J. Dhar. Stability analysisof SIR model with vaccination. J. Comp. and AppliedMath. , 17-23 (2014).[3] D. L. Chao and D. T. Dimitrov, Seasonality and the effec-tiveness of mass vaccination, Math. Biosci. Eng. ,249259 (2016).[4] H. S. Rodrigues, Application of SIR epidemiologicalmodel: new trends, arXiv:1611.02565 (2016).[5] W. H. Hamer, Epidemic disease in England – the evi-dence of variability and of persistency of type. Lancet ii,733-739 (1906).[6] R. Ross, Report on the Prevention of Malaria in Mauri-tius. London: Waterlow and Sons (1908).[7] R. Ross, An application of the theory of probabilities tothe study of a priori pathometry. – Part I. Proc. R. Soc.Lond. A , 204230 (1916).[8] R. Ross and H. P. Hudson. An application of the theoryof probabilities to the study of a priori pathometry.–PartIII. Proc. Roy. Soc. A , 225-240 (1917).[9] W. O. Kermack and A. G. McKendrick. A contributionto the mathematical theory of epidemics. Proc. Roy. Soc.A , 700-721 (1927).[10] D. G. Kendall, Discussion of Measles periodicity andcommunity size by M. S. Bartlett, J. Roy. Stat. Soc. A , 6476 (1957).[11] M. S. Bartlett, Measles Periodicity and Community Size,J. Royal Stat. Soc. A , No. 1, 48-70 (1957).[12] M. S. Bartlett, Deterministic and Stochastic Models forRecurrent Epidemics, Berkeley Symp. on Math. Statist. and Prob., Proc. Third Berkeley Symp. on Math. Statist.and Prob., Vol. 4, 81-109 (Univ. of Calif. Press, 1956).[13] W. D. Flanders and D. G. Kleinbaum, Basic Models forDisease Occurrence in Epidemiology, Int. J. Epidemiol-ogy , Issue 1, 17 (1995).[14] J. E. Amaro, The D model for deaths by COVID-19,arXiv:2003.13747v1 (2020).[15] D. S. Hui, E. Azhar, T. A. Madani et al. ., The continu-ing 2019-nCoV epidemic threat of novel coronaviruses toglobal health The latest 2019 novel coronavirus outbreakin Wuhan, China. Int. J. Infect. Dis. et al. , Inferring COVID-19 spreading ratesand potential change points for case number forecasts,https://arxiv.org/abs/2004.01105v2 (2020).[22] https://covid19.isciii.es/[23] https://dashboard.covid19.data.gouv.fr/ ∼ amaro/coronavirus/[26] https://covid19.healthdata.org/[27] Special Issue on Challenges in Modelling Infectious Dis-ease Dynamics, Edited by J. Lloyd-Smith, D. Mollison,J. Metcalf, P. Klepac and H. Heesterbeek, Epidemics ,1-108 (2015).[28] H. Heesterbeek et al. , Modeling infectious disease dynam-ics in the complex landscape of global health, Science , Issue 6227, aaa4339 (2015).[29] F. Brauer, Mathematical epidemiology: past, present,and future. Infectious Disease Modelling , 113-127(2017).[30] S. G. Krantz, P. Polyakov, A. S. R. Srinivasa Rao, Trueepidemic growth construction through harmonic analysis,J. Theor. Biology494