key: cord-0435744-b7hzmv7n
authors: Carletti, Timoteo; Fanelli, Duccio; Piazza, Francesco
title: COVID-19: The unreasonable effectiveness of simple models
date: 2020-05-22
journal: nan
DOI: nan
sha: 4669e31f2a3d8d02cc29993177ba0b7ea3d474e1
doc_id: 435744
cord_uid: b7hzmv7n

When the novel coronavirus disease SARS-CoV2 (COVID-19) was officially declared a pandemic by the WHO in March 2020, the scientific community had already braced up in the effort of making sense of the fast-growing wealth of data gathered by national authorities all over the world. However, despite the diversity of novel theoretical approaches and the comprehensiveness of many widely established models, the official figures that recount the course of the outbreak still sketch a largely elusive and intimidating picture. Here we show unambiguously that the dynamics of the COVID-19 outbreak belongs to the simple universality class of the SIR model and extensions thereof. Our analysis naturally leads us to establish that there exists a fundamental limitation to any theoretical approach, namely the unpredictable non-stationarity of the testing frames behind the reported figures. However, we show how such bias can be quantified self-consistently and employed to mine useful and accurate information from the data. In particular, we describe how the time evolution of the reporting rates controls the occurrence of the apparent epidemic peak, which typically follows the true one in countries that were not vigorous enough in their testing at the onset of the outbreak. The importance of testing early and resolutely appears as a natural corollary of our analysis, as countries that tested massively at the start clearly had their true peak earlier and less deaths overall.

In December 2019 the novel coronavirus disease SARS-CoV2 (COVID-19) emerged in Wuhan, China. Despite the drastic containment measures implemented by the Chinese government, the disease quickly spread all over the world, officially reaching pandemic level in March 2020. 1 From the beginning of the outbreak, the international community braced up in the commendable effort to rationalise the massive load of data, 2,3 painstakingly collected by local authorities throughout the world. However, epidemic spreading is a complex process, whose details depend to a large extent on the behaviour of the virus 4,5 and of its current principal vectors (we, the humans). Additionally, recorded data necessarily suffer from many kinds of bias, so that much information is intrinsically missing in the data released officially.

For example, the total number of reported infected bears the imprint of the growing size of the sampled population. This major source of non-stationarity should be properly gauged if one aims at explaining the observed trends and elaborate forecasts for the epidemics evolution. Similarly, the time series of recovered individuals that are released officially are heavily unreliable, due to the uneven ability of different countries to trace asymptomatic or mildly symptomatic patients.

In this letter, we demonstrate that overlooking these effects inevitably leads to fundamentally wrong predictions, irrespective of whether mean-field, compartment-based models 6, 7 or more sophisticate individual-based, stochastic spatial epidemic models 8, 9 are employed.

More precisely, it is important to realise that the populations of infected reported by na-tional agencies do not provide a reliable picture to correctly identify the epidemics peaks.

The primary cause resides in the marked non-stationary character of the testing activity.

However, while sources of bias due to differences in sampling and reporting policies among different countries are generally acknowledged, notably by studies applying Bayesian inference 10, 11 and causal models, 12 surprisingly the critical issue of non-stationarity appears to have been significantly overlooked in the literature on COVID-19. Nonetheless, as we show in the following, save few creditable exceptions, the reporting rate has generally increased substantially as the outbreak was gathering momentum, along with the stepping up of the number of tests conducted. Let us stress this point once more -any modelling efforts will be pointless unless non-stationarity of reporting is properly accounted for.

While infected and recovered bear the marks of the convolution with unpredictably nonstationary reading frames, it can be argued that the time series of deceased individuals provide a more reliable picture of the true, bias-free underlying dynamics of the disease. In fact, while there might have been inaccuracies in attributing deaths to COVID-19, it can be surmised that these should lead to systematic and reasonably stationary errors.

Contrary to intuition, and to an excellent degree of accuracy, a death-only analysis of the data shows that COVID-19 outbreaks in all countries fall in the simple SIR universality class, mean-field compartmental models originated by the the seminal 1927 paper by Kermack and McKendrick. 13 Importantly, this fact can only be brought to the fore in the reduced space of deaths and their time derivatives, as this allows one to filter out virtually all spurious effects associated with non-stationary reading frames, which have been systematically underestimated in the vast majority of reports.

Taken together, our analysis suggests that simple mean-field models can be meaningfully invoked to gather a quantitative and robust picture of the COVID-19 epidemic spreading.

Additionally, as a natural byproduct of our reasoning, we were able to calculate the timedependent reporting fractions for infected and recovered individuals. These provide a clear a-posteriori reading frame for the relative position of the true epidemic peak and that ob-served from the reported data, which is found to be typically shifted away artificially in the future. Furthermore, casting the SIRD dynamics as a function of the deaths time series and its derivative leads us to derive simple formulae that prove invaluable for real-time monitoring of the disease. As an important corollary of our analysis, doubts can be cast on all procedures targeted at estimating the reproduction index of the epidemics based on the raw time series of infected (and recovered), i.e. not corrected for non-stationary testing.

Within a simple mean-field compartmental model, people are divided into different species. 14 For example, in the susceptible (S), infected (I), recovered (R), dead (D) scheme (SIRD), 15 any individual in the fraction of the overall population that will eventually get sick belongs to one of the aforementioned classes. The latter represent the different stages of the disease and their size is used to model its time course. Let S 0 = S(0) be the size of the initial pool of susceptible individuals. Then the SIRD dynamics is governed by the following set of coupled

whereẊ stand for the time derivative of X = S, I, R, D. The parameter r represents the infection rate, i.e. the probability per unit time that a susceptible individual contracts the disease when she enters in contact with an infected person. The parameters a and d denote the recovery and death rates, respectively. Note that, S(t) + I(t) + R(t) + D(t) = S 0 + I 0 is an invariant of the SIRD dynamics, i.e. it remains constant as time elapses 1 In general, 1 It is reasonable to assume that R(0) = 0 as this coronavirus has never been in contact with humans before. Furthermore, S(0) I(0) O(1), which justifies the way Eqs.

(1) are written down (see also Methods).

one could attempt to calibrate model (1) against the available time series for the infected, recovered and deaths, so as to determine the best-fit estimates of r, a, d and S 0 . However, it may be argued that this procedure amounts to over-fitting, due to the intrinsic high correlation between r and S 0 in the model. More commonly, a face-value figure for S 0 is assumed from the start, e.g. the size of a country or of some specific region. However, this choice is somehow arbitrary, as the true basin of the epidemics can in principle be determined meaningfully only a posteriori.

More importantly, as argued above, there is a more fundamental limitation to any attempt of adjusting a model to the raw figures beyond specific, model-dependent issues. The cause is the time-varying bias in the reported data. Essentially, this translates into an a-priori unknown discrepancy between the true size of infected and recovered and the figures released officially. The key observation at this point is that the time series of deaths, D(t), can be regarded as the only data set that returns a faithful representation of the epidemics course.

More precisely, it appears reasonable to surmise that all unknown sources of errors on deaths can be considered as roughly stationary over the time span considered.

It is well-known that the SIR dynamics can be recast in a particularly appealing form in the space (Ṙ, R), R denoting in this case the overall class of recovered and deceased. 16, 17 When one distinguishes between the two populations D and R (the recovered), as in the SIRD dynamics (1), a similar manipulation leads to (see Methods)

where γ = dS 0 , ξ = r/(dS 0 ) and β = a + d. The differential equation (2) is equivalent to the full set of equations (1) and the associated conservation law. Overall, it provides an expedient and virtually bias-free means of testing whether the time evolution of an outbreak is governed by a SIRD-like dynamics. It appears remarkable that a process that is inherently complex and unpredictably manyfaceted can be accurately described in terms of an effective, highly simplified mean-field scheme, where many factors are deliberately omitted, such as space and different sorts of population stratification. Of course, the effectiveness of SIR schemes and simple modifications thereof are well-known. 19, 20 For example, a recent work shows that a modified SIR model where the infection rate is let vary with time can be employed to infer change points in the epidemic spread to gauge the effectiveness of specific confinement measures in Germany. 18 In an earlier work, two of the present authors had considered a similar approach to elaborate predictions of the outbreak in Italy that proved fundamentally right when compared with the evolution observed afterwards. 14 What Fig. 1 critically adds to the picture is a bias-free confirmation that simple mean-field schemes can be used meaningfully to draw quantitative conclusions.

The best-fit values of the reproduction number R 0 extracted from the above analysis should be considered as effective estimates of the overall severity of the outbreak in different countries. In particular, it can be argued that these should incorporate all sources of nonstationarity left besides the time-dependent bias associated with the evolution of the testing rate. Most likely, these amount to different degrees of (time-dependent) reduction of the infection rate r, following the containment measures. Interestingly, we find that, rather generally, the evolution of such non-stationary SIRD dynamics falls into the same universality class in the plane (Ḋ, D) as signified by f (D), Eq. (2) (see also section 3 in the supplementary material). More precisely, it turns out that a reduction of r by a factor f < 1 over a certain time window simply corresponds to a proportional rescaling of the reproduction number to a lower value, R 0 → R 0 f . A similar conclusion holds if the non-stationarity causes the opposite effect, that is, f > 1. This might be the case, for example, of overlooked hotbeds or large undetected gatherings that cause the infection to accelerate at a certain point in time.

Remarkably, whatever the direction, it can be proved that this rescaling is independent of the typical time scale over which the transition r → f r unfolds, as well as the point in time that marks its onset (see supplementary material).

An important corollary of our analysis is that the observed epidemic peak displayed by the time series of the measured infected, I M (t), typically occurs later than the true peak, the more so the more the testing rate is ramped up during the outbreak. The correct occurrence of the peak and a comprehensive understanding of how the testing activity controls the apparent evolution of the epidemics can be obtained with the following simple argument.

In general, one can define two time-dependent reporting fractions, gauging the unknown disparity between the true and the measured populations

where R M (t) refers to the reported population of recovered persons. Taking into account the first definition in Eqs. (3), differentiation of the last equation in (1) with respect to time

The apparent epidemic peak occurs whenİ M (t) = 0, whereas we know from the fact that the epidemics is governed by a SIRD-like dynamics that the true peak occurs whenD = 0, i.e.

when the number of deaths/day reaches a maximum. Thus, from Eq. (4) one immediately sees that whether the apparent peak is observed before or after the true peak depends on the sign of the rate of reporting,α I (t). More precisely, if the testing activity is steadily ramping up (α I > 0), the true peak will occur earlier than the apparent one. This is because the reported infected will have a maximum forD < 0, i.e. past the maximum of theḊ.

Conversely, if the testing rate decreases, this will anticipate the apparent peak, giving the false impression that the worst might be over while the true number of infected is in fact still increasing. We find that this analysis applies to all countries considered in this paper (see also supplementary material), whereby either the former or the latter scenarios are invariably Figure 2 : The derivative of the reporting fraction of infected controls the apparent epidemic peak. Top panels: fraction of true recovered reported, α R (t) (grey shaded regions) and fraction of true infected reported, α I (t), computed as described in the text. The shaded regions correspond to an assumed mortality m = 0.012 ± 0.005. 21 The thick lines mark averages on the last 15 % portions of the data ( α I,R T ), and first 40 % ( α I 0 ), while the dashed lines denote the overall average values of α I,R computed over the whole available time span. Bottom panels. The reported infected are compared to the rescaled deaths/day,

the infected that would have been reported if the reporting fraction had remained constant at its initial value α I (0). Technically, δ = d/α I (0) is obtained as the slope δ of a linear fit ofḊ vs I M in the very early stage of the outbreak. The shaded regions correspond to the least-square errors found on δ. As soon as the testing rate starts varying appreciably, the two trends I 0 (t) and I M (t) diverge. If the testing rate is ramped up, the infected reported seem to flag an increase in the outbreak severity (Italy). Conversely, if the testing rate slows down (Germany), the reported infected underestimate the true peak. In particular, in this case the peak can be alarmingly anticipated. China shows a marginal behaviour, where the two peaks coincide, consistently with a remarkably stable reporting fraction of infected. observed. (3), it is not difficult to obtain

where β is estimated from the fits of the SIRD model in the plane (Ḋ, D). Direct inspection of Fig. 2 confirms the above reasoning on the relative position in time of the apparent and true peaks. For example, it can be seen that tests have increased exponentially in Italy, starting from the end of March, which caused the apparent peak to occur nearly a month after the true one. Conversely, the testing activity has been vigorous in Germany at the onset of the outbreak, but appears to have slowed down substantially as the epidemics progressed. As a consequence, the apparent peak has been reported while the epidemics was still growing. In both cases, our estimate is consistent with the fact that only 15 % of infected persons are effectively reported. China seems to provide an instance of the marginal caseα I = 0. Despite an early ramping up, the reporting fraction of infected appears to be remarkably stable, which is consistent with the co-occurence of the true and the apparent peaks. According to the data, the testing activity has been resolute since the beginning (only about half of the infected eluding tracing), and has culminated with a scenario compatible with virtually all infected individuals being traced.

The calculations reported in Fig. 2 also confirm that counting instances of recovery is a much subtler task. Typically, countries seem to have needed some adjustment time to set up robust counting protocols, as it is apparent from the large fluctuations that are often observed in the trends of α R at early times (see also supplementary material). Another interesting feature that emerges clearly from the comparison of α I and α R is the average recovery time, which can be inferred for example by direct inspection of the testing activity in Italy. The curve α R (t) appears to be shifted in the future by about 10 days with respect to α I (t), which is consistent with the present consensus estimate of an average time for recovery of two weeks. Figure 3 : The true epidemic peak occurs within 30 days after the first ten deaths, the earlier the higher the initial testing rate. When shifted so that they coincide at the day where ten new deaths were first reported, theḊ time series reveal that the true epidemic peak occurs inevitably rather early in the outbreak, visibly not later than one month after the first ten casualties reported. Earlier vigorous testing is obviously key in containing the outbreak (see also discussion in the text).

Our analysis shows how to pinpoint the true peak of the epidemics -the maximum of the deaths/day time series,Ḋ. It is thus natural to investigate how the different true peaks compare among different countries once these time series have been brought to a common time origin. This is illustrated in Fig. 3 for a number of European countries. It can be clearly appreciated that the apogee of the outbreak in different countries does not occur later than about one month after the first ten deaths reported. It is interesting to observe that in countries where the testing rate was vigorous at the onset of the outbreak (large α I 0 ), e.g.

Austria, Switzerland, Belgium, Germany 2 , the peak seems to occur earlier with respect to countries that have started ramping up their testing capacity later into the outbreak, such as explanation that early tracking of infected individuals enables more effective isolation from the start, thus leading to a less severe evolution of the outbreak as there simply are less contagious people able to propagate the disease. It is also interesting to observe that peak position and height show some correlation with the initial reproduction numbers measured in the (Ḋ, D) plane (see again Fig. 1 and Tables S1 and S2 in the supplementary material).

With some notable exceptions (again Netherlands and Denmark), countries that totalised higher deaths/day at peak tend to be characterised by higher values of R 0 . Overall, effective and timely testing seems to have been the key to successful containment strategies, as it is by now widely advocated, 22 as opposed to muscular lockdowns imposed later into the outbreak.

In summary, the observed behavior of the early 2020 COVID-19 outbreak has raised many diverse and puzzling issues in different countries. Many of them have adopted drastic containment measures, which had some effect on the spread of the virus but also left large sectors of the population world-wide grappling with the daunting perspective of pernicious economic recessions. On their side, scientists from all disciplines are having a hard time identifying recurring patterns that may help design and tune more effective response strategies for the next waves. In this letter we have shown that the spread of the epidemics caused by this elusive 4 pathogen falls into the universality class of SIR models and extensions thereof.

Based on this fact, an important corollary of our analysis is that early testing is not only key in making sense of the reported figures, but also in controlling the overall severity of the outbreak in terms of casualties.

It is intrinsically difficult to disentangle the effects imposed by containment measures from the natural habits of a given community in terms of social distancing. For example,

Italy is recording three times as many deaths per inhabitant than Switzerland as we write.

However, while it is true that Switzerland seemed to have tested nearly ten times as much initially (see supplementary material), it is also true that that country has never been put in full lockdown as it was decided in Italy. At the same time, Swiss people are known for their law-abiding attitude and social habits that make (imposed or suggested) social distancing certainly easier to maintain than in southern countries such as Italy.

The key message of this letter is that simple models should not be dismissed a priori in favour of more complex and supposedly more accurate schemes, which unquestionably go into more detail but at the same time rely necessarily on a large number of parameters that it is hard to fix unambiguously. All in all, we can state with certainty that vigorous early testing and accurate and transparent self-evaluation of the testing activity, assuredly virtuous attitudes in general, should be considered all the more a priority in the fight against COVID-19.

The data used in this paper span the interval 1/22/20 -5/16/20 and have been retrieved from the github repository associated with the interactive dashboard hosted by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, Baltimore, USA. 23

In this section we derive Eq. with ξ = r/(dS 0 ). Furthermore, since R(0) = 0, one has from Eqs. (1) that R(t) = (a/d)D(t).

Taking into account that by definition I(t) =Ḋ/d, the conservation law, S(t) + I(t) + R(t) + D(t) = S 0 + I 0 , can be recast in the following form Supplementary information D = δI M . The errors reported have been computed by constructing confidence ellipses at 90 % confidence level (γ, β, ξ), and as standard least-squares errors for the linear fit to extract the slope δ. Fig. 4 . Table 2 ). A direct proportionality is observed, in agreement with relation (S.1). The vertical and horizontal bars denote the errors reported in Table 2 , calculated according to standard error propagation from the formulae used to compute the two quantities. The orange dashed line is the arbitrary proportionality relation D = S 0 /100, drawn to provide an order of magnitude of the ratio between deaths at peak and effective size of the infection basin.

In this section we derive some useful relations about key quantities at the infection peak, namely the time t at which the number of infected persons, I(t), reaches its maximum value,

. This can be determined from the conditionİ(t) = 0. From the fourth equation of (1) in the main text, we getD(t) = 0, then using the definition of f (D) we can write

where D = D(t). Being D a monotonic increasing function we have,Ḋ(t) > 0 for all t, hence the latter condition is equivalent to

where the stand for the derivative of f with respect to the variable D. This equation can be straightforwardly solved to give

where R 0 = r/(a + d) (see also Methods in the main text). We are now able to compute the remaining key quantities at the peak time t. Indeed, the relations developed in the Methods section in the main text imply

Let us observe that these quantities, as well as D, scale proportionally to the system size, S 0 , as one can appreciate from Fig. 4 where we report D as a function of S 0 . The latter is computed as S 0 = γ/d (see main text). Here d = δα I (0), where δ is the slope of a linear fiṫ D vs I M at short times and γ is a parameter of the function f (D) obtained from a fit in the plane (Ḋ, D) (see Fig. 1 in the main text).

In the main text we have shown (see Fig.1 where N = S 0 + I 0 ≈ S 0 . Then (see Eq.

(2) in the main text and Methods)

where we have used the relation R 0 = ξγ/β. The universal function is unimodal, φ(0; R 0 ) = 0

for all values of R 0 and there is a single positive root,x(R 0 ), which converges to 1 for It is easy to see that φ (x R 0 ; R 0 ) → −1 as R 0 → ∞. In order to explore this question, we consider a modified version of the SIRD model (Eqs.

(1) in the main text), where the infection rate r is let vary with time. 14 More precisely, if some containment measures were enforced at time t * , causing their effect on a typical time ∆t, we may take

where r 0 is the initial, pre-lockdown infection rate that defines the true initial reproduction Whatever the effect of the non-stationarity introduced by r(t), the effect in the rescaled plane (x , x) turns out to be extremely simple and apparently robust with respect to the choice of the characteristic time scales t * , ∆t and also of the specific functional form that The effective reproduction numbers are plotted in the right panel against the corresponding reduction of infection rate, f . These seem to be described by the same universal linear reduction, R eff 0 (f ) = R 0 f , irrespective of the characteristic time scales that govern the non-stationarity and even of the specific function the drives the reduction from r 0 to f r 0 < r 0 (see text). describes the reduction from r 0 to f r 0 . This is illustrated in Fig. 6 for a choice of f in the interval [0.1, 3]. The evolution of outbreaks generated from a given root one, that is, one with a given infection rate r(0) = r 0 , can always be described by the same one-parameter function φ(x; R eff 0 ) with an effective reproduction number R eff 0 that only depends on f . Moreover, R eff 0 for f < 1 is simply the t = 0 value reduced by the same factor, that is

Remarkably Egypt Morocco Algeria Figure 9 : A positive derivative of the reporting fraction of infected causes a delay of the apparent epidemic peak with respect to the true one -III. Top panels: fraction of true recovered reported, α R (t) (grey shaded regions) and fraction of true infected reported, α I (t), computed as described in the text. The shaded regions correspond to an assumed mortality m = 0.012 ± 0.005. 21 The thick lines mark averages on the last 15 % portions of the data ( α I,R T ), and first 40 % ( α I 0 ), while the dashed lines denote the overall average values of α I,R computed over the whole available time span. Bottom panels. The reported infected are compared to the rescaled deaths/day, I 0 (t) def =Ḋα I (0)/d, i.e. the infected that would have been reported if the reporting fraction had remained constant at its initial value α I (0). Technically, δ = d/α I (0) is obtained as the slope δ of a linear fit oḟ D vs I M in the very early stage of the outbreak. The shaded regions corresponds to the least-square errors found on δ.

Finland Canada Figure 10 : A negative derivative of the reporting fraction of infected deceptively anticipates the apparent epidemic peak with respect to the true one. Top panels: fraction of true recovered reported, α R (t) (grey shaded regions) and fraction of true infected reported, α I (t), computed as described in the text. The shaded regions correspond to an assumed mortality m = 0.012 ± 0.005. 21 The thick lines mark averages on the last 15 % portions of the data ( α I,R T ), and first 40 % ( α I 0 ), while the dashed lines denote the overall average values of α I,R computed over the whole available time span. Bottom panels. The reported infected are compared to the rescaled deaths/day, I 0 (t) def =Ḋα I (0)/d, i.e. the infected that would have been reported if the reporting fraction had remained constant at its initial value α I (0). Technically, δ = d/α I (0) is obtained as the slope δ of a linear fit oḟ D vs I M in the very early stage of the outbreak. The shaded regions corresponds to the least-square errors found on δ.

Israel Denmark Figure 11 : The outbreak in a reduced number of marginal cases is consistent with an overall flat reporting fraction of infected. Top panels: fraction of true recovered reported, α R (t) (grey shaded regions) and fraction of true infected reported, α I (t), computed as described in the text. The shaded regions correspond to an assumed mortality m = 0.012 ± 0.005. 21 The thick lines mark averages on the last 15 % portions of the data ( α I,R T ), and first 40 % ( α I 0 ), while the dashed lines denote the overall average values of α I,R computed over the whole available time span. Bottom panels. The reported infected are compared to the rescaled deaths/day, I 0 (t) def =Ḋα I (0)/d, i.e. the infected that would have been reported if the reporting fraction had remained constant at its initial value α I (0). Technically, δ = d/α I (0) is obtained as the slope δ of a linear fit ofḊ vs I M in the very early stage of the outbreak. The shaded regions corresponds to the least-square errors found on δ.

New coronavirus outbreak: Framing questions for pandemic prevention

With COVID-19, modeling takes on life and death importance

Unique epidemiological and clinical features of the emerging 2019 novel coronavirus pneumonia (COVID-19) implicate special control measures

Why do some COVID-19 patients infect many others, whereas most don't spread the virus at all? Science

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Multiscale mobility networks and the spatial spreading of infectious diseases

Why is it difficult to accurately predict the COVID-19 epidemic? Infectious Disease Modelling

Modeling the COVID-19 outbreaks and the effectiveness of the containment measures adopted across countries

The Stanford Encyclopedia of Philosophy, summer

Containing Papers of a Mathematical and Physical Character

Analysis and forecast of COVID-19 spreading in China, Italy and France

Data-based analysis, modelling and forecasting of the COVID-19 outbreak

Forecasting the outcome and estimating the epidemic model parameters from the fatality time series in COVID-19 outbreaks

Priesemann, V. Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation

How lethal is the novel coronavirus, and how many undetected cases there are? The importance of being tested

Covid-19 mass testing facilities could end the epidemic rapidly

An interactive web-based dashboard to track COVID-19 in real time