A parsimonious description and cross-country analysis of COVID-19 epidemic curve
AA parsimonious description and cross-countryanalysis of COVID-19 epidemic curves
Kristoffer Rypdal and Martin RypdalDepartment of Mathematics and StatisticsUiT – The Arctic University of Norway
Abstract
In a given country, the cumulative death toll of the first wave of the COVID-19epidemic follows a sigmoid curve as a function of time. In most cases, the curve is welldescribed by the Gompertz function, which is characterized by two essential param-eters, the initial growth rate and the decay rate as the first epidemic wave subsides.These parameters are determined by socioeconomic factors and the countermeasuresto halt the epidemic. The Gompertz model implies that the total death toll dependsexponentially, and hence very sensitively, on the ratio between these rates. The re-markably different epidemic curves for the first epidemic wave in Sweden and Norwayand many other countries are classified and discussed in this framework, and theirusefulness for the planning of mitigation strategies is discussed.
At the time of writing, early August 2020, the COVID-19 pandemic is expanding exponen-tially in many countries, particularly in the tropics and the southern hemisphere. On theother hand, in Southeast Asia, Europe, and North America, the first wave of infection ismostly over. In some countries, e.g., the United States, Israel, Iran, and Spain, a secondand stronger wave develops. For the first wave, the general shape of the epidemic curves,the daily number of confirmed cases or COVID-19 related deaths, is one of rapid growthfollowed by a slower decay. However, even though this general characteristic is ubiquitous,the total death toll per million inhabitants in comparable countries varies by more than anorder of magnitude. For instance, by August 5, 2020, Sweden had registered 569 deaths permillion, while neighboring Norway only 47.The debate about the causes of this pronounced variability between countries has involveda plethora of explanations [1] based on correlation techniques such as multiple regression[2], or on models that involve statistical modelling of the dynamics of transmission [3].Suggested predictors (explaining variables) in such analyses have encompassed demographic,socioeconomic, and health related characteristics, and strength and timing of border closures1 a r X i v : . [ q - b i o . P E ] A ug nd social distancing measures. Some success in assessing the effect of non-pharmaceuticalinterventions in 11 European countries was obtained in [3], and the effect of social distancingmeasures in Sweden, Denmark and Norway was studied by a considerably simpler approachin [4]. The three countries are similar in terms of health care, language, culture, climate,economy, and institutional framework. Moreover, due to the simultaneous import of infectedcases via ski tourists returning from Austria and Italy in early March 2020, it is plausible thatcommunity spread of the disease started at approximately the same time in these countries.On the other hand, while Norway and Denmark followed the general lockdown of most ofthe European Union countries on March 12, 2020, Sweden followed a different path, largelyusing voluntary measures [5].On a world scale, cross-country comparative analysis is extremely challenging for at leasttwo reasons. One is that the relevant predictors vary across the world and may attain anunmanageable number. Another is that the same problem could be pertain to the predictands(explained variable). How do we best describe the impacts of the disease in a way thatallows comparison between countries in different parts of the world? Some studies employregression analysis to a sample containing the majority of the world’s countries, but considerthe accumulated COVID-19 cases and/or deaths up to a certain date as predictand [6]. Thisvariable may be very misleading, however, since the virus was not introduced simultaneouslyin all countries. In principle, we need to predict and compare the full epidemic curvesacross countries, but in order to extract information from from such comparison, we needa simplified characterization of these curves, i.e., we need to reduce the information in thecurves to a few numbers. This paper is about how we can obtain and make use of such acharacterization.Our main finding is that the epidemic curves for COVID-19 related deaths for mostcountries with a reliable reporting system are surprisingly well described by the so-calledGompertz growth model [7]. This model contains only two essential free parameters to bedetermined by the data. These parameters determine the initial relative growth rate of thenumber of daily deaths and the decay rate of this number as the first epidemic wave comesto an end. The mathematical properties of the Gompertz model imply that the total deathtoll in the first wave is given by the exponential of the ratio between these two rates. Thisratio can be considered as a measure of the skewness of the curve for daily deaths. Theexponential dependence implies that the death toll is strongly dependent on this skewnessparameter. Countries with rapid initial growth and slow later decay suffer the higher deathtoll, but there are variations among countries with a similar number of deaths. For instance,while Spain experienced very fast initial growth followed by a rapid decay due to a veryforceful lockdown, Swedish death numbers rose and decayed considerably slower. Othercountries, like Norway, experienced slower initial growth than Sweden, and faster decay, andthereby obtained death numbers one order of magnitude lower. Since the initial growth ratein Norway was so much lower than in Spain, Norway could obtain these results with slowerdecay, i.e., with a considerably milder lockdown, than Spain.In this paper, we make a detailed comparison of epidemic curves for Sweden and Norway,including computation and plotting the respective time evolution of the relative growth2ates and time-dependent reproduction numbers. We further fit the Gompertz model to theepidemic curves (deaths) for a large number of countries that have completed the first wave,and plot the model parameters in a diagram highlighting their differences. The two ratesare predictands that correlate in different ways to various predictors of the epidemic curve,and via these predictands they determine the total death toll in the first wave. Regressionsthat determine the coefficients of these predictors will be useful in the design of optimalmitigation strategies for upcoming waves of the epidemic. Data employed in this paper are downloaded from Our World in Data [8]. Their data forCOVID-19 related deaths are retrieved from the European Centre for Disease Preventionand Control. As explained by the COVID-19 Health System Response Monitor [9], figuresmay vary among countries and may complicate cross-country analysis. This problem is mostserious for the headline figures presenting the most recent day-to-day data. We use datapublished at the end of July, but only data up to the first week of July. In Section 4 we willdiscuss the implications of such uncertainties for the conclusions of the analysis.The official figures from China suffer from 40% discontinuous increase in death numberson April 17. The official explanation is that home deaths before that date were included. Inour analysis, we account for this by adjusting the numbers before April 17 by a factor 1.4,thus removing the discontinuity.
The Gompertz model is a special case of the class of logistic growth models [10]. This generalclass are solutions J ( t ) to the nonlinear, separable differential equation dJdt = γ ( J ) J, (1)where t is the time coordinate and J is a quantity undergoing monotonic growth that satu-rates at the carrying capacity J ∞ in the limit t → ∞ . In the Gompertz model the instanta-neous relative growth rate γ ( J ) = d t J/J = d t (ln J ) is (see Appendix A), γ ( J ) = γ ∞ ln (cid:18) J ∞ J (cid:19) , (2)and Eq. (1) can be integrated to yieldln J ( t ) = ln J ∞ − ln (cid:18) J ∞ J (cid:19) exp ( − γ ∞ t ) , (3)3r equivalently, J ( t ) = J ∞ (cid:18) J J ∞ (cid:19) exp ( − γ ∞ t ) . (4)Note that in this limit γ ( J ) diverges as J →
0, so Eq. (1) should be understood as an initialvalue problem with J (0) = J and an initial growth rate γ = γ ∞ ln (cid:18) J ∞ J (cid:19) . (5)The logarithmic form of Eq. (3) shows that the Gompertz model describes a growth wherethe logarithm of J converges exponentially towards the limit ln J ∞ , i.e., the difference J ∞ − J decays with the decay rate γ ∞ . The growth curve is completely determined by the initialvalue J , the carrying capacity J ∞ , and the asymptotic decay rate γ ∞ . These features makeit natural to plot J ( t ) in a logarithmic plot.From Eq. (3) we have ln( J ∞ /J ) = exp ( − γ ∞ t ) ln( J ∞ /J ), and Eq. (2) then yields γ as afunction of time, γ ( t ) = γ exp ( − γ ∞ t ) . (6)This is an alternative way of expressing equations (1) and (2), which shows that the relativegrowth rate as a function of time decays exponentially with rate γ ∞ for all t > γ ( t ) depends only on γ and γ ∞ . Eq. (5) shows that γ is determined by the three model parameters J , J ∞ , and γ ∞ , and hence can replace anyof these as a model parameter. The parameters J ∞ , and γ ∞ are more fundamental than theothers, however, since the latter depend on the choice of the time origin. Eq. (4) impliesthat a translation of the time variable, t → t − t , corresponding to a shift of the time origin,leads to a shift of the initial value; J → J ∞ ( J /J ∞ ) exp ( − γ ∞ t ) , and Eq. (6) leads to a shiftof the initial growth rate; γ → γ exp ( γ ∞ t ). A natural choice is to choose the time originto be the first day the observed value of J exceeds one death per million inhabitants.In the fitting procedure we do not fix J = 1, thus allowing the modeled value J tobe slightly different from one and from the observed value at t = 0. Eq. (5) can now berewritten to give us the total death toll in terms of the rates γ and γ ∞ , J ∞ = J exp (cid:18) γ γ ∞ (cid:19) , (7)where J is close to one. The significance of Eq. (7) is that the accumulated death tollover the first epidemic wave depends exponentially on the ratio of the growth rate γ at thetime when the number of deaths exceeds one per million, to the asymptotic decay rate γ ∞ as the epidemic burns out. A more asymmetric epidemic curve; rapid rise and slow decay,intuitively leads to a larger number of deaths, but the exponential dependence on the ratioof rates suggests that the sensitivity of the death toll with respect to countermeasures issurprisingly high.The monotonically increasing sigmoid shape of J ( t ) suggest that it models the evolutionof accumulated quantities like the cumulative number of infected individuals or cumulative4
20 40 60 80 10051050100 time ( days ) C u m u l a t i v edea t h s ( days ) D ea t h s pe r da y (a) (b) Γ ! γ " Figure 1: (a): Evolution of cumulative deaths per million J ( t ) in a log-plot expressed as aGompertz function. (b): The evolution of the daily death number d t J ( t ) in a log-plot withthe tangents at t = 0 and t = ∞ marked. The slopes of these tangents are Γ and γ ∞ ,respectively.deaths. In the literature one also frequently deals with the daily numbers, i.e. with thederivative d t J . Fig. 1(a) shows an example of J ( t ) given by the Gompertz function, andFig. 1(b) its derivative, both in logarithmic plots. In Fig. 1(b)Γ( t ) = d t ln( d t J ) = γ ∞ (cid:18) e − γ ∞ t ln (cid:18) J ∞ J (cid:19) − (cid:19) = γ ∞ (cid:20) ln (cid:18) J ∞ J ( t ) (cid:19) − (cid:21) . (8)is the relative growth rate of the daily number d t J which is a skew, bell-shaped func-tion. Γ( t ) is a monotonically decreasing function starting out at the positive value Γ = γ ∞ (ln ( J ∞ /J ) − d t J -curve, and converging towards thenegative value − γ ∞ as t → ∞ . From Eq. (8) we find an alternative to Eq. (7) J ∞ = J exp (cid:18) γ ∞ (cid:19) = ( J e ) exp (cid:18) Γ γ ∞ (cid:19) , (9)and it follows that Γ = γ − γ ∞ . (10)Thus, we observe that the total death toll depends exponentially on the ratio Γ /γ ∞ , whichis the absolute value of the ratio of the initial and the final slopes of the curve in Fig. 1(b).Hence this ratio represents a measure of the skewness of this curve.The parameters J ∞ and γ ∞ are intrinsic in the sense that they do not depend on the timeorigin, while γ and Γ do exhibit such a dependence. This is reflected by the inconvenientfact that the estimated J is different for each country and is determined by the choice wehave made for the time origin (the first day the number of deaths per million exceeds one).For comparison of countries, it would be more correct to use the growth rate γ at the time t when J ( t ) = 1. This growth rate is found by putting J = 1 in Eqs. (2) and (8), yielding γ = γ ∞ ln J ∞ , Γ = γ − γ ∞ , (11)5uch that Eqs. (7) and (9) reduce to, J ∞ = exp (cid:18) γ γ ∞ (cid:19) = exp (cid:18) γ ∞ (cid:19) = e exp (cid:18) Γ γ ∞ (cid:19) , (12)or, by taking the logarithm, ln J ∞ = γ γ ∞ = 1 + Γ γ ∞ (13) R ( t ) In Appendix A the Gompertz model is derived heuristically from the SIR dynamical modelwith some auxiliary assumptions. In the SIR model, the variable J ( t ) is the cumulativenumber of infections in a country’s population. This quantity cannot be observed directly,because the majority of infected individuals are never tested for the virus. The number ofcases confirmed by tests is not a reliable proxy, because testing regimes change over thecourse of the epidemic, and vary among countries as well. The number of COVID-19 relateddeaths is also unreliable as a proxy, since infection death rates vary from one country toanother, but as a measure of the time development of the cumulative number of infectionsit is probably the best statistic that is generally available. In this paper, we are interestedin the overall shape of the epidemic curve, and not the time for epidemic onset or detailsdepending on the distribution of the time between infection and death. For most countriesthe shape of the curves for infections and deaths are very similar in the sense that the formercan be obtained by a trivial time translation and re-scaling of the latter J ( t ) → s d J ( t + t d ) , (14)where the scaling factor s d in Scandinavia is of the order 150 and the time lag between infec-tion and death t d is about 20 days. After performing this re-scaling of the curve of observeddeaths, we interpret the result as the curve of cumulative infection cases that appears in theSIR model (see Appendix A). In that model we encounter the time-dependent reproductionnumber R ( t ), which is more commonly used in epidemiology than the instantaneous relativegrowth rate γ ( t ). Consider now the the linearized SIR equations (30) and (31) for the cumu-lative number of infected, J ( t ), and the number of active transmitters of the infection, I ( t ),as described in Appendix A. Introducing the reproduction number R ( t ) = β ( t ) /α , theseequations take the form dJdt = α R I, (15) dIdt = α ( R − I. (16)Here, the reproduction number is typically reduced with time due to countermeasures andthe gradual removal of superspreaders from the susceptible population, while the recovery6ate α remains constant. By solving Eq. (16) and inserting into Eq. (15), the latter can bewritten, dJdt = αI R exp (cid:32) (cid:90) t α ( R − dt (cid:48) (cid:33) . (17)Taking the logarithm of this equation and differentiating with respect to t , yields the differ-ential equation d R dt + α R − α R − Γ( t ) R = 0 , (18)where Γ( t ) is given by Eq. (8). The equation can be solved numerically with the properinitial condition. In the asymptotic limit t → ∞ , Γ( t ) → − γ ∞ and R ( t ) → R ∞ , where R ∞ = 1 − γ ∞ α . (19)From Eq. (17) we find that for small t the growth rate is γ ( t ) ≈ α ( R ( t ) − t = 0, we have, by using Eq. (2) and the convention J = 1, R = 1 + γ α = 1 + γ ∞ α ln J ∞ . (20)Thus, in terms of R and R ∞ , Eq. (7) becomes J ∞ = exp (cid:18) R − − R ∞ (cid:19) . (21) The essential methodological idea in this work is to fit a mathematical model for the epidemiccurve to an observed time series; in this case the cumulative number of COVID-19 relateddeaths. This curve has a sigmoid shape with near-exponential growth in the initial growthphase. In such situations a usual least-square fitting procedure will give a poor representationof the initial growth because large relative errors will give small contributions to the absoluteleast-square error. Much better representation is obtained by fitting the logarithm of themodel to the logarithm of the data. This is what we do here, by fitting the function for ln J ( t )given by Eq. (3) to the logarithm of the cumulative death per million time series startingat the first day this number exceeds one. We use the built-in fitting routine FindFit in Mathematica J ∞ , γ ∞ , and J .We have typically used about 100 days of the time series for this fitting, but have adjustedthis manually for each country to make sure that we only include the first wave of theepidemic. We have only included countries which have clearly completed a first wave, sothat we can reliably estimate the decay rate γ ∞ . For this reason, we have not analyzedmany countries in the tropics or the southern hemisphere. The data for the 73 countries wehave analyzed, and the fits to them, are shown in Figs. 9-14. The fits are surprisingly good,and demonstrate the usefulness of the Gompertz model to give an analytical representationof these data. 7 ��� (cid:1) �� Apr May Jun Jul Aug51050100500 date c u m u l a t i v edea t h s pe r m illi on Apr May Jun Jul0.000.050.100.150.200.250.300.35 date γ ( t ) (a) (b) Figure 2: (a): Red bullets represent cumulative COVID-19 related deaths per millioninhabitants in Norway from March 23 and onwards, and the full curve is the Gompertzcurve J ( t ) fitted to these data. The blue bullets and curve are the corresponding for Sweden(shifted 3 days forwards). Note that the plot is logarithmic, and that cumulated death tollper million in Sweden in early July is 12 times that in Norway. (b): The relative growthrate γ ( t ) = d t (ln J ) as given by Eq. (6) for Norway (red) and Sweden (blue). These growthrates are the slopes of the curves in (a). In this section we first present a detailed comparison of Sweden and Norway in order to high-light some of the useful information that can be drawn from the analytical representation ofthe epidemic curves through the fitted Gompertz functions. Next we estimate the Gompertzmodel parameters from 73 countries and classify them in terms of total death toll and ratiobetween rates of initial growth and later decay.
Fig. 2(a) shows the development of cumulative deaths per million inhabitants in Norwayand Sweden in a logarithmic plot, and the corresponding fitted Gompertz function. Swedendata converge to a total death toll which is more the ten times higher than that of Norway.The relative growth rates γ ( t ) = d t ln J ( t ) given as the slope of each curve are plotted inFig. 2(b). Eq. (6) shows that these growth rates decay exponentially from an initial value γ towards zero at a rate γ ∞ . Observe that the gap between Swedish and Norwegian deathnumbers continues to grow throughout April and May (see also the growth of the ratio ofthese numbers in Fig. 4(b)), consistent with the higher Swedish growth rate throughout thisperiod.Fig. 3(a) shows the derivative d t J ( t ) for Sweden an Norway in logarithmic plot. Thegrowth rate Γ( t ) of this derivative is the slope of these curves, and is plotted in fig. 3(b).It is given by Eq. (8) and decays exponentially from Γ = γ − γ ∞ to − γ ∞ . As shown byEq. (12) the total death toll J ∞ is determined by the ratio of the initial and final slope, i.e.,the ratio Γ /γ ∞ which characterizes the skewness of the curve. The steeper initial growthin Sweden signified through the higher Γ obviously contributes to the total death toll, but8 ��� (cid:1) �� Apr May Jun Jul Aug0.51510 date da il y dea t h s pe r m illi on Apr May Jun Jul - Γ ( t ) (a) (b) Figure 3: (a): Evolution of daily deaths per million, dJ/dt , computed numerically from Eq.(12) for Norway (red) and for Sweden (blue). (b): The relative growth rate Γ( t ) = d t (ln d t J )for Norway (red) and for Sweden (blue).so does also also the lower rate on the decaying slope in June and July.Before we try to quantify these contributions, let us look at the reconstructed evolutionof the reproduction number R ( t ) for the two countries by solving Eq. (18) with Γ( t ) shifted17 days backwards in time to account for the average delay between infection and death, and α − = 5 days representing the average infectious period. The initial value R is given byEq. (20). The 17-days delay is uncertain, so the shape of R -curves shown in Fig. 4(a) maybe more accurate than the absolute dating. The important observation is that the relativelysmall, but consistent, difference between the reproduction number in Sweden and Norwaythroughout the entire epidemic wave has been sufficient to produce a tenfold higher numberof deaths in Sweden.Fig. 4(b) shows that the ratio of the cumulative deaths in Sweden and Norway increasesstrongly throughout the entire epidemic wave, which indicates that the higher R in Swedenduring the decaying phase has had a strong effect on the total death toll, i.e., people havecontinued dying in Sweden in May and June, long after the death rate was close to zero inNorway.Further insight into the differences of the epidemic curves between countries can beobtained by the graphics demonstrated in Fig. 5. Here we have marked the point (Γ , γ ∞ )for Sweden and Norway and represented the function ln J ∞ = 1 + Γ /γ ∞ in a density plot.The diagonal lines represent isolines for constant J ∞ . We observe that Sweden is located onthe iso-line where J ∞ = 520 and Norway on the line where J ∞ = 49. In addition, it is seen,as already noticed, that Sweden exhibits a considerably higher initial growth rate Γ thanNorway, and also a lower decay rate γ ∞ . The main value of the diagram, however, is that itallows us to explore the effect on the death toll of hypothetical action in the two countriesthat did not materialize.Judging from the final death tolls, the Swedish path does not seem as a very attractiveoption. But what if the countries had chosen to follow the Swedish path initially, with stronggrowth in the early phase and then followed the Norwegian path with fast decay in the latephase? The end state would then be the purple cross which is located on the isoline with9 ��� (cid:1) �� Apr May Jun0.00.51.01.52.02.53.0 date R ( t ) Apr May Jun Jul0246810 date J S w eden / J N o r w a y (a) (b) Figure 4: (a): Evolution of R ( t ) computed numerically from Eq. (12) for Norway (red)and for Sweden (blue). (b): Evolution of ratio between accumulated deaths in Sweden andNorway.totally 224 deaths per million. Another option would be to follow the Norwegian path withslow initial growth followed by the Swedish path with slow decay. This would lead to theblue cross which lies on the 82 deaths per million line. This suggests that reduction of theinitial growth is the most effective way to bring the total death toll in the first wave down.More precisely, the lower initial growth in Norway reduced the number of deaths by a factorsix compared to Sweden, and the faster decay in the late phase by nearly another factor two.It should be mentioned that the hypothetical paths discussed above may not be realizablein practice since there may be dynamical constraints that do not allow Γ and γ ∞ to bevaried independently of each other. However, as described in Section 3.2 by studying a largenumber of countries we find a wide range of possible states in the (Γ , γ ∞ )-plane, and nosign of “forbidden” areas within this range. Figs. (6) shows the same as Fig. (5), where the data for Sweden and Norway has beensupplemented by 71 other countries. The epidemic curves and their Gompertz fits are shownin Fig. (9) – (14). Fig. (6) shows that the death toll J ∞ in the first wave covers a range ofmore than two orders of magnitude, from just above 3 in China to more than 800 in Belgium.Ordered by increasing death toll (see Fig. 7(c) for the full ordered list), we have identified afew countries in the legend. China, New Zealand, Slovakia, South Korea and Japan belongto the group of countries with less than 10 deaths per million. In the group with ∼ deaths, we have Iceland, Norway, Finland, Austria, Denmark, and Germany. And in groupwith ∼ deaths, United States, Sweden, Spain, Italy, United Kingdom, and Belgium.From Fig. 7(c) one could be tempted to conclude that the epidemic curves are verysimilar for countries with similar death toll, but Fig. 6 and Figs. 7(a,b) show that this isnot so. The magnitudes of the growth rates vary considerably among countries with similar J ∞ = exp (Γ /γ ∞ ). Some countries within one group have a rapid initial rise Γ which is10 ! = ! ! = ! ! = ! ! = Norway Sweden
Figure 5: Density plot of J ∞ (Γ , γ ∞ ) with some isolines marked, along with the points forSweden and Norway. The crosses mark the end states for paths where the inital and finalphase of the two countries are mixed. 11 ! = ! ! = ! ! = (cid:1)(cid:2)(cid:1) (cid:1)(cid:2)(cid:3) (cid:1)(cid:2)(cid:4) (cid:1)(cid:2)(cid:5) (cid:1)(cid:2)(cid:6)(cid:1)(cid:2)(cid:1)(cid:4)(cid:1)(cid:2)(cid:1)(cid:6)(cid:1)(cid:2)(cid:1)(cid:7)(cid:1)(cid:2)(cid:1)(cid:8)(cid:1)(cid:2)(cid:3)(cid:1) (cid:1) (cid:1) (cid:1) (cid:1) AustriaBelgiumBrazilChinaDenmarkFinlandGermanyIcelandItalyJapan New ZealandNorwaySlovakiaSouth KoreaSpainSwedenUnited KingdomUnited States
Figure 6: Density plot of J ∞ (Γ , γ ∞ ) with some isolines marked, along with the points for73 countries. The legend shows the positions for some selected countries.compensated by a rapid fall (high rate γ ∞ ) in the late phase, while others have a slow initialgrowth and a slower decay, resulting in approximately the same death toll. As mentionedin the introduction, Sweden experienced considerably slower initial growth than Spain, butalso achieved much slower decay than Spain due to weak measures to contain the epidemiccompared to Spain’s radical lock-down. The result is that the two countries are located veryclose to each other in Fig. 7(c), but are strongly separated along the same J ∞ isoline in Fig.6. Mathematically, Fig. 6 locates each country in a 3D plot which specifies the three Gom-pertz parameters (Γ , Γ ∞ , ln J ∞ ). Hence, it contains all the information contained in theGompertz description of the epidemic curve. Fig. 7 can be conceived as the projections ofthe points in Fig. 6 onto the three axes. The main purpose of this paper has been to extract objective information from the epidemiccurves for COVID-19 related deaths on the country level by utilizing the observational factthat the curves for the first wave are quite accurately represented by the Gompertz functionfor many countries of interest. The results may be subject to different interpretations and12 a)(b)(c) B u r k i na F a s o N i ge r J a m a i c a G eo r g i a C h i na M a l a ys i a S i ngapo r e M a li A u s t r a li a S ou t h K o r ea T un i s i a S o m a li a B ang l ade s h C a m e r oon A f ghan i s t an P a k i s t an P h ili pp i ne s G abon M o r o cc o J apan U k r a i ne N e w Z ea l and B u l ga r i a S l o v a k i a C uba B e l a r u s M on t eneg r o La t v i a R u ss i aL i t huan i a G r ee c e M e x i c o C y p r u s B o s n i aand H e r z ego v i na A l ban i a P o l and M o l do v a C r oa t i a B r a z il P e r u D o m i n i c an R epub li c K u w a i t P ana m a E c uado r S e r b i a H unga r y P ue r t o R i c o R o m an i a F i n l and I c e l and N o r w a y I s r ae l C z e c h R epub li c I r an T u r k e y C anada S l o v en i a G e r m an y E s t on i a D en m a r k P o r t uga l A u s t r i a S w eden U n i t ed S t a t e s S w i t z e r l andLu x e m bou r g I r e l and U n i t ed K i ngdo m F r an c e It a l y N e t he r l and s B e l g i u m S pa i n (cid:1) ( per day ) N e w Z ea l and S l o v a k i a C h i na E s t on i a T un i s i a A u s t r i a C z e c h R epub li c A l ban i a S pa i n S l o v en i a M o r o cc o C uba A u s t r a li a I s r ae l D en m a r k N e t he r l and s B e l g i u m I c e l and B u r k i na F a s oLu x e m bou r g S w i t z e r l and P o r t uga l G e r m an y M a l a ys i a C y p r u s I r e l and T u r k e y S e r b i a F r an c e It a l yJ apan P ue r t o R i c o G r ee c e N o r w a y U n i t ed K i ngdo m N i ge r I r an U n i t ed S t a t e s G eo r g i a C r oa t i aLa t v i a F i n l and S w eden S o m a li a S ou t h K o r ea H unga r y R o m an i a J a m a i c a M on t eneg r o C anada P o l andL i t huan i a S i ngapo r e D o m i n i c an R epub li c P ana m a B o s n i aand H e r z ego v i na P h ili pp i ne s K u w a i t E c uado r M o l do v a M a li B e l a r u s B u l ga r i a R u ss i a U k r a i ne P e r u C a m e r oon B r a z il G abon M e x i c o P a k i s t an B ang l ade s h A f ghan i s t an (cid:1) (cid:1) ( per day ) B u r k i na F a s o N i ge r C h i na J a m a i c a M a l a ys i a G eo r g i a T un i s i a A u s t r a li a N e w Z ea l and S i ngapo r e S l o v a k i a S ou t h K o r ea M o r o cc o S o m a li a C uba J apan M a li A l ban i a P h ili pp i ne s C y p r u s La t v i a G r ee c e M on t eneg r o C a m e r oonL i t huan i a C r oa t i a S e r b i a C z e c h R epub li c I c e l and I s r ae l P o l and U k r a i ne B u l ga r i a P ue r t o R i c o G abon N o r w a y E s t on i a S l o v en i a B o s n i aand H e r z ego v i na B e l a r u s T u r k e y D o m i n i c an R epub li c H unga r y F i n l and B ang l ade s h A u s t r i a R o m an i a I r an D en m a r k P a k i s t an G e r m an y R u ss i a P ana m a K u w a i t P o r t uga l M o l do v aLu x e m bou r g S w i t z e r l and A f ghan i s t an C anada N e t he r l and s U n i t ed S t a t e s E c uado r I r e l and F r an c e S w eden P e r u It a l y B r a z il S pa i n M e x i c o U n i t ed K i ngdo m B e l g i u m (cid:1) / (cid:2) (cid:1) Figure 7: (a): Estimated Γ for the 73 countries ranked from lower to higher values. (b):Estimated γ ∞ ranked from high to low values. (c): Estimated Γ / Γ ∞ = ln J ∞ ranked fromlower to higher values. 13aise many questions that we cannot answer in this paper. In this section, we will limitourselves to present some subjective views on why the Gompertz model seems to fit thesedata so well and some viewpoints on the information presented in Figs. 6 and 7. By taking the time derivative of Eq. (3), the Gompertz model can be written in the form γ ( t ) = d t JJ = d t ln J ( t ) = γ ∞ ln (cid:18) J ∞ J (cid:19) exp( − γ ∞ t ) = γ exp( − γ ∞ t ) , (22)which is identical to Eq. (6). The relative growth rate γ represents how much one unitof J grows per unit time. It is closely related to the reproduction number R , which canbe interpreted as the average number of new infections caused by one infected individual.If the disease transmission mechanisms remain constant over time, including the density ofsusceptible individuals, γ will remain constant, and the growth of J ( t ) will be exponential. Ifa large fraction of the population becomes infected, and eventually immune, herd immunitywill appear as a nonlinearity in the SIR-equations (see Appendix A). This effect has not yetappeared on the country level in the COVID-19 epidemic, so changes in the transmissionmechanism must cause a reduction of γ ( t ) over time. One effect that is seen in all epidemicseven without societal action is that individuals who are particularly active in spreadingthe disease (so-called super-spreaders) catch the disease early and are removed from thesusceptible population. As time goes, this reduces the effective reproduction number and γ . In addition, society acts in complex ways to resist the disease. If the total effect of thisresistance on the rate of change of γ is proportional γ itself (which is a common property ofa dissipative force), the equation for γ would take the form dγdt = − ηγ, (23)which is what we will find by taking the time derivative of Eq. (22) with the resistance η = γ ∞ . Thus, we have arrived at a straightforward interpretation of the Gompertz modelas a mathematical description of an unstable system where a quantity J ( t ) naturally growsexponentially in the absence of dissipation. The dissipative force in this system does notact on J ( t ) itself, but rather on its growth rate. Society does not take action in response tothe death toll itself but to its growth. A physics analog could be the force on an electricallycharged body that emits electromagnetic radiation. This force resists the acceleration, therate of change of the velocity, unlike an ordinary friction force, which resists the velocity. The interpretation of Figs. 6 and 7 is sensitive to how the sample of countries has beenselected. A complete picture will emerge as more countries complete the first wave. We hadexcluded countries that are not well beyond the first peak in daily deaths, or countries that14ad entered a second wave much before the first was completed. Countries that have notreached the first peak are Argentina, India, South Africa, and Colombia, and countries thathave entered an early second wave are Iran, Saudi Arabia, Iraq, and Chile.We have also excluded countries that we strongly suspect have unreliable or irregularreporting or epidemic curves that, by visual inspection, are not well fitted by a sigmoidfunction. However, we emphasize that no country has been excluded from the sample basedon the results of the analysis.The relative effect of the two parameters on the total death toll is seen by traversing thesample in Fig. 6 in the vertical and horizontal direction. We observe that the distribution inthe γ ∞ -direction is much wider for small Γ than for large Γ . In the vertical direction, J ∞ varies two orders of magnitude among countries with Γ < .
15, i.e., there is a high variabilityof the decay rate γ ∞ among countries with low initial growth. Some countries, like China,New Zealand, Slovakia, South Korea, and Japan, have taken strong action throughout thefirst wave, despite low initial growth of the death numbers. On the other hand, othercountries like Brazil, have used relatively low initial growth as an excuse for non-action,resulting in a very high death toll. For countries with Γ > .
15 the range of γ ∞ becomesnarrower with increasing Γ . In this group, we find most Western-European countries andthe United States. Finland, Norway, Iceland, Germany, Denmark, and Austria have Γ inthe lower end of this range and J ∞ ∼ or less, while for Sweden, United States, UnitedKingdom, Italy, Belgium, and Spain, Γ is higher and J ∞ ∼ . The overall impression isthat the variability of the death toll in Western Europe and the United States is largely dueto variations of the initial growth rate, and to a lesser extent due to variations of the laterdecay. There are exceptions, however. For instance, Sweden and Austria have almost thesame Γ , but γ ∞ on opposite tails of the distribution, yielding about ten times higher deathtoll in Sweden.In the rest of the world, the picture is more mixed. Some low-income countries, withlow initial growth rates, still end up with a high death toll because the decay rate is alsolow. On the other hand, high-income countries in South-East Asia suffer few deaths due toa combination of low initial growth combined with moderate or high decay rate.The exceptionally low Γ for China is due to the strong confinement of the epidemic tothe Hubei province, which constitutes only 4.3 percent of the Chinese population. If wehad treated Hubei as a country, the death toll per million would have been 23 times higher,and Γ would have increased by a factor ln 23 ≈ is defined as the growth rate at the time the death toll exceeds one permillion, and this time comes earlier if the population is considered to be only that of Hubei.This observation underscores that geographic isolation of the epidemic to limited regionswithin a country is crucial in reducing the initial growth rate and the total death toll. Thesuccess of the Chinese strategy in limiting the first wave is the effective isolation of Hubeifrom the rest of China and the very strict lockdown within the province.A caveat of this entire discussion is that the death rates reported from the various coun-tries may be inaccurate. Systematic under-reporting will influence the estimate of Γ , but notthat of γ ∞ . At present, we have not been able to make systematic corrections to these figures15or all countries. Corrections based on figures for excess mortality is a possibility for somecountries, but such figures do not exist for many of countries for which the official COVID-19 death rates are least reliable. Excess mortality may also give rise to under-estimationof COVID-19 deaths in countries with effective interventions, because these interventionsreduce mortality from other diseases. The huge variability of the death toll among countries in the first wave of the COVID-19epidemic is puzzling. Many causal explanations have been suggested in the media, but on aworld scale, little systematic work has been published in the peer-reviewed literature. Thepresent work does not attempt a causal explanation, but seeks a quantitative characterizationof the epidemic curve in terms of as few parameters as possible. We demonstrate that thisis possible in terms of the Gompertz function, which describes a three-parameter sigmoidcurve. One of these parameters is redundant because it only determines the time origin thatdetermines the point in time when the death toll per million inhabitants exceeds one. Thetwo parameters that determine the total death toll is the initial relative growth rate of thenumber of daily deaths and the decay rate of this number as the first wave subsides. Thesalient property of the Gompertz function is that the total death toll of the first wave dependsexponentially on the ratio of those rates. It hence describes the exceptional sensitivity ofthe death toll to the values of these rates.We illustrate these features by applying the analysis to compare the effect of initial growthrates and later decay rates of the neighbouring countries Norway and Sweden, where morethan ten times higher death toll in Sweden is attributed to a higher early death rate and, toa lesser extent, a lower decay rate in the subsiding phase. It is also shown to be related toa somewhat higher reproduction number in Sweden throughout the entire wave.The Gompertz function is shown to be well fitted to the death data for 73 countriesaround the world, and plots of the Gompertz parameters for these countries show thatthey are broadly distributed in this parameter space, suggesting that there are many routesto a given death toll. A task for future work is to find out more about commonalitiesbetween countries that are close to each other in the parameter space, and about featuresthat distinguish countries that appear similar but end up at different locations in this space.The most straightforward approach would be some sort of multiple regression analysis bywhich the growth rate Γ and the decay rate γ ∞ are linked to socioeconomic variables andintervention measures.At the time of writing, many countries are entering the second epidemic wave. It will beinteresting to investigate if this wave can be analyzed by the same method and to study ifparameters change from the first to the second wave. A more complete data set, involvingmore countries, will also become available, facilitating a more complete classification of theepidemic curves of the COVID-19 pandemic.16 ppendixA From the SIR model to the Gompertz model Let N be the total population, S the number of susceptible individuals, and I the numberof infected. We then have, dSdt = − β ISN , (24) dIdt = β ISN − αI, (25)where β is the rate by which the infection is being transmitted, and α is the rate by whichthe infected are isolated from the susceptible population. Decoupled from these, we have theequation for the number R of individuals that are quarantined, deceased, and recovered; dRdt = αI. (26)These three equations constitute the SIR model. A.1 Nonlinear evolution of an “old” disease
During the development of an epidemic, health authorities typically report the accumulatednumber of registered infected J = N − S . The equations for J and I then take the form dJdt = β (cid:18) − JN (cid:19) I, (27) dIdt = (cid:18) β − α − β JN (cid:19) I. (28)When one studies the flare-up of an old pathogen, there may be a considerable herd immunity,i.e., J may be of the same order of magnitude as N . Then, the peak growth rate of I isattained when J = J s = N (cid:18) − αβ (cid:19) , (29)and the saturation and decay of the epidemic is a nonlinear effect. A.2 Evolution of an epidemic due to a ”new” pathogen
If the pathogen is new, there is no immunity in the population from the start, and if weassume that only a small fraction of the population is affected by the epidemic, we have J (cid:28) N for all times and we can linearize; dJdt = βI, (30)17 Idt = ( β − α ) I. (31)Initial growth, and later decay, of I now requires that α and β evolve in time, and that β − α changes sign from positive to negative at a certain time t s . Since J ( t ) is monotonicallyincreasing with t , there is a bijective mapping between t and J , and we can think of α , β ,and I as functions of J rather than t . A.3 The reproductive number
If we define R ( J ) = β/α as the reproductive number, Eqs. (24) and (25) reduce to dIdJ = 1 − R ( J ) . (32)It is not unreasonable to assume that R ( J ) is a monotonically decreasing function, sincemeasures to reduce the transmission rate will lower β as J increases and measures to diagnoseand quarantine infected individuals will increase α . I ( J ) has its maximum at J = J s = J ( t s )for which R ( J s ) = 1. This is also the inflection point for J ( t ); (cid:18) d Jdt (cid:19) t s = 0 . (33) A.4 Reduction to a generalized logistic growth model
Since I = I ( J ), the equation dJdt = − βI, (34)can be written as dJdt = γ ( J ) J, (35)where γ ( J ) = β IJ , (36)is the relative growth rate for J . Eq. (34) has the general form of a nonlinear growth modelwith a growth rate depending on the growing variable itself. A.5 Richards’ growth model J ( t ) takes some sort of sigmoid form, which means that γ ( J ) is monotonically decreasingfrom γ ( J = 0) = γ (0) to γ ( J ∞ ) = 0, where J ∞ = J ( t → ∞ ). A common form is to express γ ( J ) by a function depending on the three parameters γ (0) , J ∞ , and ν ; γ ( J ) = γ (0) (cid:20) − (cid:18) JJ ∞ (cid:19) ν (cid:21) , (37)18 ��� (cid:1) �� J / J ∞ γ / γ ( ) t J / J ∞ (a) (b) n =1/8 n =1/4 n =1/2 n =1 n =2 n =4 n =8 n =8 n =1 n =1/2 n =1/4 n =1/8 Figure 8: (a): Normalized growth rate γ/γ (0) as a function of
J/J ∞ for varying ν . For ν = 1we have the symmetric logistic growth curve. Increasing ν > J s closer to the limiting value J ∞ . Decreasing ν < J . (b): NormalizedRichards growth curve J/J ∞ as a function of time t for varying ν .and hence Eq. (34,) reduces to Richards’ equation, dJdt = γ (0) (cid:20) − (cid:18) JJ ∞ (cid:19) ν (cid:21) J. (38)Note that γ (0) = γ ( J = 0) is the growth rate for J = 0, not the growth rate at time t = 0.The latter is γ = γ ( t = 0) = γ (0) (cid:18) − (cid:18) J J ∞ (cid:19) ν (cid:19) , (39)The normalized growth rate γ ( J/J ∞ ) /γ (0) is plotted for varying ν in Fig. (8)(a), and thecorresponding Richards curves in Fig. (8)(b). The solution satisfying the initial condition J ( t = 0) = J is J ( t ) = J ∞ (cid:0) J ∞ /J ) ν − e − γ (0) νt (cid:1) /ν , (40)leaving us with a model with four free parameters to be determined by fitting the model tothe time series for the death toll; γ (0) , J ∞ , ν , and J . In the asymptotic limit t → ∞ , J ( t )converges monotonically towards the limit J ∞ , and the difference from this limit takes theform J ∞ − J ( t ) → ν (cid:20)(cid:18) J ∞ J (cid:19) ν − (cid:21) e − γ (0) νt , (41)hence it decays exponentially, and so does the derivative J (cid:48) ( t ) (the daily number of deaths),at a rate γ ∞ = γ (0) ν. (42)19 .6 The Gompertz limit; ν → , γ (0) ν → γ ∞ < ∞ It is observed from Figure 8 that the limit ν → γ ( t ), followed by a slower decay towards zero as J → J ∞ . By fittingthe Richards model to time series for countries that have gone through the first wave of theepidemic, we typically find small values for the exponent ν . This suggests that for countrieswith a reasonably rapid response to the increasing death tolls it may be reasonable to lookat this limit, and hence reduce the number of free parameters by one. The resulting model isthe one suggested by Gompertz [7]. Using the relation lim ν → ν − [1 − ( J/J ∞ ) ν ] = ln ( J ∞ /J ),it is easily shown that in this limit Eq. (37) reduces to γ ( J ) = γ ∞ ln (cid:18) J ∞ J (cid:19) , (43)and Eq. (40) to J ( t ) = J ∞ (cid:18) J J ∞ (cid:19) exp ( − γ ∞ t ) , (44)Note that in the Gompertz limit, γ ( J ) diverges as J →
0, while for finite ν (Richards’model) it converges towards γ (0) . The important feature, however, is not the growth duringthe zoonotic stage, or the stage when the epidemic is dominated by imported cases, butrather the stages when the pathogen is transmitted in the population. By defining the timeorigin t = 0 as the time when J exceeds a certain threshold J , and considering only t ≥ γ = γ ( J ) = γ ∞ ln (cid:18) J ∞ J (cid:19) , (45)which is completely determined by the model parameters; J , J ∞ , and γ ∞ . B Gompertz fits to 73 countries
Figs. 9-14 presents the Gompertz fits to the 73 countries with parameters plotted in Figs. 6and 7. 20igure 9: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 21igure 10: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 22igure 11: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 23igure 12: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 24igure 13: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 25igure 14: Shows cumulative COVID19 deaths (black) and the estimated Gompertz curvesfor different countries (red). 26 eferences [1] Sung, W.W.Y., and Kaplan, R.M. Why do countries’ COVID-19 death ratesvary so much? Multiple factors at play: testing capacity, case defi-nitions, age distribution, preparedness. MEDPAGE TODAY May 15 , , doi: 10.1016/j.eclinm.2020.100464[3] Flaxman S., Mishra S., Gandy, A. et al. Estimating the effects of non-pharmaceuticalinterventions on COVID-19 in Europe. Nature, June 8 .[4] Juranek, S., and Zoutman, F.T. The Effect of Social Distancing Measures on IntensiveCare Occupancy: Evidence on COVID-19 in Scandinavia. Discussion Papers ,Norwegian School of Economics, Department of Business and Management Science.[5] Nature, April 30 , Vol. 580 , d41586-020-01098-x.pdf[6] Farzanegan, R.F., Feizi, M., and Gholipour, H.F. Globalization and outbreak ofCOVID-19: an empirical analysis. CESifo working paper no. 8315 May .[7] Gompertz, B. On the nature of the function expressive of the law of human mortal-ity, and on a new mode of determining the value of life contingencies. PhilosophicalTransactions of the Royal Society of London , , 513585.[8] Our World in Data https://ourworldindata.org/coronavirus-data[9] Karanikolos, M., and McKee M. How Comparable is COVID-19 Mortal-ity Across Countries? COVID-19 Health System Monitor, June 4 .https://analysis.covid19healthsystem.org/index.php/2020/06/04/how-comparable-is-covid-19-mortality-across-countries/[10] Tsoularis, A. Analysis of logistic growth models, Res. Lett. Inf. Math. Sci. , ,23–46.[11] Verhulst, P.F. Notice sur la loi que la population suit dans son accroissement, Corr.Math. Physics ,10