Inferring change points in the COVID-19 spreading reveals the effectiveness of interventions

As COVID-19 is rapidly spreading across the globe, short-term modeling forecasts provide time-critical information for decisions on containment and mitigation strategies. A main challenge for short-term forecasts is the assessment of key epidemiological parameters and how they change when first interventions show an effect. By combining an established epidemiological model with Bayesian inference, we analyze the time dependence of the effective growth rate of new infections. Focusing on the COVID-19 spread in Germany, we detect change points in the effective growth rate that correlate well with the times of publicly announced interventions. Thereby, we can quantify the effect of interventions, and we can incorporate the corresponding change points into forecasts of future scenarios and case numbers. Our code is freely available and can be readily adapted to any country or region.

Introduction. When faced with the outbreak of a novel epidemic like COVID-19, rapid response measures are required by individuals as well as by society as a whole to mitigate the spread of the virus. During this initial, time-critical period, neither the central epidemiological parameters, nor the effectiveness of interventions like cancellation of public events, school closings, and social distancing are known.
Rationale. As one of the key epidemiological parameters, we infer the spreading rate λ from confirmed COVID-19 case numbers at the example of Germany by combining Bayesian inference based on Markov-Chain Monte-Carlo sampling with a class of SIR (Susceptible-Infected-Recovered) compartmental models from epidemiology. Our analysis characterizes the temporal change of the spreading rate and, importantly, allows us to identify potential change points and to provide short-term forecast scenarios based on various degrees of social distancing. A detailed description is provided in the accompanying paper, and the models, inference, and predictions are available on github. While we apply it to Germany, our approach can be readily adapted to other countries or regions.
Results. In Germany, interventions to contain the outbreak were implemented in three steps over three weeks: Around March 9, large public events like soccer matches were cancelled. On March 16, schools and childcare facilities as well as many non-essential stores were closed. One week later, on March 23, a far-reaching contact ban ("Kontaktsperre"), which included the prohibition of even small public gatherings as well as the further closing of restaurants and non-essential stores, was imposed by the government authorities.
From the observed case numbers of COVID-19, we can quantify the impact of these measures on the disease spread (Fig. 0). Based on our analysis, which includes data until April 21, we have evidence of three change points: the first changed the spreading rate from λ 0 = 0. With this third change point, λ transitioned below the critical value where the spreading rate λ balances the recovery rate µ, i.e. the effective growth rate λ * = λ − µ ≈ 0 ( Fig. 0, gray traces). Importantly, λ * = 0 presents the watershed between exponential growth or decay. Given the delay of approximately two weeks between an intervention and first inference of the induced changes in λ * , future interventions such as lifting restrictions warrant careful consideration.
Our detailed analysis shows that, in the current phase, reliable short-and long-term forecasts are very difficult as they critically hinge on how the epidemiological parameters change in response to interventions: In Fig. 0 already the three example scenarios quickly diverge from each other, and consequently span a considerable range of future case numbers. Thus, any uncertainty on the magnitude of our social distancing in the past two weeks can have a major impact on the case numbers in the next two weeks. Beyond two weeks, the case numbers depend on our future behavior, for which we have to make explicit assumptions. In the main paper we illustrate how the precise magnitude and timing of potential change points impact the forecast of case numbers (Fig. 2).
Conclusions. We developed a Bayesian framework to infer central epidemiological parameters and the timing and magnitude of intervention effects. Thereby, the efficiency of political and individual intervention measures for social distancing and containment can be assessed in a timely manner. We find evidence for a successive decrease of the spreading rate in Germany around March 6 and around March 15, which significantly reduced the magnitude of exponential growth, but was not sufficient to turn growth into decay. Our analysis also shows that a further decrease of the spreading rate occurred around March 23, turning exponential growth into decay. Future interventions and lifting of restrictions can be modeled as additional change points, enabling short-term forecasts for case numbers. In general, our analysis code may help to infer the efficiency of measures taken in other countries and inform policy makers about tightening, loosening and selecting appropriate rules for containment. until April 1 -two c.p. until April 9 -three c.p.

Forecasts with limited data (95% CI)
Data, reported cases

Start of interventions:
* ** *** Inferring change points in the COVID-19 spreading reveals the effectiveness of interventions INTRODUCTION During the initial outbreak of an epidemic, reliable short-term forecasts are key to estimate required medical capacities, and to inform and advice the public and decision makers [1]. During this initial phase, three tasks are of particular importance to provide time-critical information for crisis mitigation: (1) establishing central epidemiological parameters, such as the basic reproduction number, that can be used for short-term forecasting; (2) simulating the effects of different possible interventions aimed at the mitigation of the outbreak; (3) estimating the actual effects of the measures taken -to rapidly adjust them and to adapt short-term forecasts. Tackling these tasks is challenging due to the large statistical and systematic errors that are present during the initial stages of an epidemic with its low case numbers. This is further complicated by the fact that mitigation measures are taken rapidly, while the outbreak unfolds, but they take an effect only after an a priori unknown delay. To obtain reasonable parameter estimates for short-term forecasting and policy evaluation despite these complications, any prior knowledge available needs to be integrated into modeling efforts to reduce uncertainties. This includes knowledge about basic mechanisms of disease transmission, recovery, as well as preliminary estimates of epidemiological parameters from other countries, or from closely related pathogens. The integration of prior knowledge, the quantitative assessment of the remaining uncertainties about epidemiological parameters, and the principled propagation of these uncertainties into forecasts is the domain of Bayesian modeling and inference [2,3].
We combine the SIR model (and generalizations thereof) with Bayesian parameter inference and augment the model by a time-dependent spreading rate. The time dependence is implemented via potential change points that reflect changes in the spreading rate driven by governmental interventions. Based on three distinct measures taken in Germany, we detect three corresponding change points from reported COVID-19 case numbers. Already on April 1 we had reported evidence for the first two change points, and predicted the third one [30]. Now, with data until April 21, we have evidence for all three change points. First, the spreading rate decreased from 0. 43 [2,9]). This matches the cancellation of large public events such as trade fairs and soccer matches. Second, the spreading rate decreased further to 0.15 (CI [0.12, 0.20]) initiated around March 15 (CI [13,17]). This matches the closing of schools, childcare facilities, and non-essential stores. Third, the spreading rate decreased further to 0.09 (CI [0.06, 0.13]) initiated around March 23 (CI [20,25]). This corresponds well to the strict contact ban, which was announced on March 22. While the first two change points were not sufficient to switch from growth of novel cases to a decline, the third change point brought this crucial reversal.
Our framework is designed to infer the effectiveness of past measures and to explore potential future scenarios along with propagating the respective uncertainties. In the following, we demonstrate the potential impact of timing and magnitude of change points, and report our inference about the three past governmental interventions in Germany. Our framework can be readily adapted to any other country or region. The code (already including data sources from many other countries), as well as the figures are all available on Github [31].

BACKGROUND: INFERENCE OF CENTRAL EPIDEMIOLOGICAL PARAMETERS AND THE EFFECTS OF INTERVENTIONS
In order to simulate the general effect of different possible interventions on the spread of COVID-19 in Germany, we first focus on the initial phase of the outbreak when no serious mitigation measures were implemented. In the absence of interventions, an epidemic outbreak can be described by SIR models with constant spreading rate (Methods). In Germany, first serious interventions occurred around March 9 and affected the case reports with an observation delay (a combination of incubation period with median 5-6 days [32]) and a test delay (time until doctor is visited plus test-evaluation time) that we assume to be both about 2-3 days. Hence, in order to infer central epidemiological parameters, we consider as the initial phase the time period from March 2 to March 15. In order to simulate the effect of different possible interventions, we then model the effects of interventions as change points in the spreading rate (Methods).
Bayesian inference for central epidemiological parameters during the initial phase of the outbreak We perform Bayesian inference for the central epidemiological parameters of an SIR model using Markov-Chain Monte Carlo (MCMC) sampling (Fig. 1). The central parameters are the spreading rate λ, the recovery rate µ, the reporting delay D and the number of initially infected people I 0 . We chose informative priors based on available knowledge for λ, µ and D, and we chose uninformative priors for the remaining parameters (Methods). Also, we intentionally kept the informative priors as broad as possible such that the data would constrain the parameters (Fig. 1).
As median estimates, we obtain for the spreading rate λ = 0.41, µ = 0.12, D = 8.6, and I 0 = 19 (see Fig. 1 C-H for the posterior distributions and the 95% credible intervals). Converted to the basic reproduction number R 0 = λ/µ, we find a median R 0 = 3.4 (CI [2.4, 4.7]), which is consistent with previous reports that find median values between 2.3 and 3.3 [18,33,34]. Overall, the model shows good agreement with both new cases ( Fig. 1 A) and cumulative cases ( Fig. 1 B) that show the expected exponential growth (linear in log-lin plot). The observed data are clearly informative about λ, I 0 and σ (indicated by the difference between the priors (gray line) and posteriors (histograms) in Fig. 1 D,E,F). However, µ and D are largely determined by our prior choice of parameters (histograms match gray line in Fig. 1 C,H). This is to be expected for the initial phase of an epidemic outbreak, which is dominated by exponential growth.
In order to quantify the impact of possible interventions, we concentrate on the effective growth of active infections before and after the intervention. As long as the number of infections and recoveries are small compared to the population size, the number of active infections per day can be approximated by an exponential growth (Fig. 1A,B) with effective growth rate λ * = λ − µ (see Methods). As a consequence, λ and µ cannot be estimated independently. This is further supported by a systematic scan of the model's log-likelihood in the λ-µ space that shows an equipotential line for the maximum likelihood ( Fig. 1 I). This strongly suggests that the growth rate λ * is the relevant free parameter with a median λ * = 28% (Fig. 1 G). The control parameter of the dynamics in the exponential phase is thus the (effective) growth rate: If the growth rate is larger than zero (λ > µ), case numbers grow exponentially; if the growth rate is smaller than zero (λ < µ), the recovery dominates and the new cases decrease. The two different dynamics (supercritical and subcritical, respectively) are separated by a critical point at λ * = 0 (λ = µ) [35]. We simulate different, hypothetical interventions following the initial phase in order to show that both, the amount of change in behavior (leading to a change in spreading rate λ, Fig. 2 A) and the exact timing of the change (Fig. 2 B) determine the future development. Hypothetical interventions build on the inferred parameters from the initial phase ( Fig. 1, in particular median λ 0 = 0.41 and median µ = 0.12) and were implemented as change points in the spreading rate from the inferred λ 0 to a new value λ 1 . With such a change point, we model three potential scenarios of public behavior: (I) No social distancing; Public behavior is unaltered and the spread continues with the inferred rate (λ 1 = λ 0 with median λ 1 = 0.41 > µ). (II) Mild social distancing; The spreading rate decreases to 50%, (λ 1 = λ 0 / 2 with median λ 1 = 0.21 > µ). Although people effectively reduce the number of contacts by a factor of two in this second scenario, the total number of reported cases continues to grow alongside scenario (I) for the time period of the reporting delay D (median D = 8.6 from initial phase, see below for a more constrained estimation). Also, we still observe an exponential increase of new infections after the intervention becomes effective, because the growth rate remains positive, λ * 1 = λ 1 − µ > 0. (III) Strong social distancing; Here, the spreading rate decreases to 10%, (λ 1 = λ 0 / 10 with median λ 1 = 0.04 < µ). The assumptions here are that contacts are severely limited, but even when people stay at home as much as possible, some contacts are still unavoidable. Even under such drastic policy changes, no effect is visible until the reporting delay D is over. Thereafter, a quick decrease in daily new infections manifests within two weeks (delay plus change point duration), and the total number of cases reaches a stable plateau.  Only in this last scenario a plateau is reached, because here the growth rate becomes negative, λ * < 0, which leads to decreasing numbers of new infections.
Furthermore, the timing of an intervention matters: Apart from the strength of an intervention, its onset time has great impact on the total case number ( Fig. 2 B,C). For example, focusing on the strong intervention (III) -where a stable plateau is reached -the effect of advancing or delaying the change point by just five days leads to more than a three-fold difference in cumulative cases.
While we find that the timing of an intervention has great effect on case numbers, the duration over which the change takes place has only minor effect -if the intervals of change are centered around the same date. In Fig. 2 C we illustrate the adjustment of λ 0 → λ 1 for mild social distancing with durations of 14, 7 and 1 day(s). The change point duration is a simple way to incorporate variability in individual behavior, and is not linked to the reporting delay D. As an interesting effect, a sudden change in the spreading rate can lead to a temporary decrease of new case numbers, despite the fact that the effective growth rate remains positive at all times.

RESULTS
In order to model real-world data, we further refined the SIR model. Most importantly, we account for systematic variations of case reports throughout the week, in particular lower case numbers on the weekend, by explicitly modeling a weekly reporting modulation (see Methods). Indeed, model comparisons confirm that models with this correction outperform those without (see Table S2). In the supplemental material, we further generalize our model to include an explicit incubation period (SEIR-like, Fig. S3) that yields results consistent with our main model.
We incorporate the effect of governmental interventions into our models by introducing flexible change points in the spreading rate (see Methods). During the COVID-19 outbreak in Germany, governmental interventions occurred in three stages from (i) the cancellation of large events with more than 1000 participants (around March 9), through (ii) closing of schools, childcare facilities and the majority of stores (in effect March 16), to (iii) the contact ban and closing of all non-essential stores (in effect March 23). The aim of all these interventions was to reduce the (effective) growth rate λ * = λ − µ. As soon as the growth rate becomes negative (λ * < 0), the number of new confirmed infections decreases after the respective reporting delay.   Table S2).
Detecting change points in the spreading rate -and quantifying the amount of change as quickly as possiblebecomes a central modeling challenge when short-term forecasts are required. To address this challenge, we assume an initial spreading rate λ 0 (the exponential growth phase, cf. Fig. 1) and up to three potential change points motivated by the German governmental interventions: The first change point (λ 0 → λ 1 ) is expected around March 9 (t 1 ) as a result of the official recommendations to cancel large events. A second change point (λ 1 → λ 2 ) is expected around March 16 (t 2 ), when schools and many stores were closed. A third change point (λ 2 → λ 3 ) is expected around March 23 (t 3 ), when all non-essential stores were closed, and a contact ban was enacted. We model the behavioral changes that are introduced at these change points to unfold over a few days ∆t i , but the changes in duration can be partly compensated by changes in the onset time t i (see Fig. 2 C, scenarios). We chose priors for all parameters based on the information available to us up to March 28 (see Methods). In addition, we performed a sensitivity analysis by employing wider priors in the supplemental material (Figs. S5, S6, S7, Table S2), which yielded consistent results. On March 28, the data were already informative about the first change point, and thereby helped to inform our forecast scenarios.
The inferred parameters for the models with change points are consistent with the inferred parameters from the exponential onset phase (Figs. 1, 3 & Figs. S1, S2). In particular, all estimated λ 0 -values from models with multiple change points are compatible with the value of the model without change points (during the exponential onset phase, λ 0 = 0.41, CI [0.32, 0.51], assuming a stationary λ until March 15, Fig. 1 E). Also the scale factor σ and the number of initial infections I 0 for the models with change points are consistent with the initial model inference during the exponential onset phase (Fig. 1 D,F).
The models with two or three change points fit the observed data better than those with fewer change points The models with three change points describe the data better than models with fewer change points, as indicated by the leave-one-out (LOO) cross-validation-based Bayesian model comparison [36] (lowest LOO-score in Table I). However, the LOO-scores of the model with two and three change points differ by less than one standard error. This originates from an extended duration of the second change point in the two-change-point model, which partially captures the effect of the third intervention. As expected, the models with none or a single change point have LOO-scores that are at least one standard deviation higher (worse) than those of the best models, and we will not consider them further.
When comparing our inference based on three change points to the number of confirmed cases, we find them to largely match (Fig. 3 B,C). The dominant periodic change in the daily new reported cases (Fig. 3 B) is well described by the weekday modulation. In addition to the periodic change, the daily new case numbers also reflect the fairly sudden change of the spreading rate at the change points (cf. Figs. 2 and S4 for the effect of change points without the modulation). Most importantly, the cumulative effect of change points manifests in an overarching decay in new case numbers that is visible after April 5 and follows the third change point (with reporting delay).

Change point detection quantifies the effect of governmental interventions on the outbreak of COVID-19 in Germany
Ideally, detected changes can be related to specific mitigation measures, so that one gains insights into the effectiveness of different measures (Fig. 3). Indeed, we found clear evidence for three change points in the posterior distributions of the model parameters: First, λ(t) decreased from λ 0 = 0.43 (with 95% credible interval, CI [0.35, 0.51]) to λ 1 = 0.25 (CI [0.20, 0.30]). The date of the change point was inferred to be March 6 (CI [2,9])]; this inferred date matches the timing of the first governmental intervention including cancellations of large events, as well as increased awareness. After this first intervention, the (effective) growth rate λ * (t) = λ(t) − µ decreased by more than a factor 2, from median λ 0 − µ = 0.3 to median λ 1 − µ = 0.12, given that the recovery rate was inferred as µ = 0.13 (CI [0.09, 0.18]). Second, λ(t) decreased from λ 1 = 0.25 to λ 2 = 0.15 (CI [0.12, 0.20]), which is larger than our prior assumption. The date of the change point was inferred to be March 15 (CI [13,17])]; this inferred date matches the timing of the second governmental intervention including closing schools and some stores. After this second intervention, the median growth rate became λ * (t) = λ 2 − µ = 0.02 ≈ 0 and is thus in the vicinity of the critical point, yet still positive. The first two interventions in Germany thereby mitigated the spread by drastically reducing the growth rate, but the spread of COVID-19 remained exponential. Third, λ(t) decreased from λ 2 = 0.15 to λ 3 = 0.09 (CI [0.06, 0.13]). The date of the change point was inferred to be March 23 (CI [20,25])]; this inferred date matches the timing of the third governmental intervention including contact ban and closing of all non-essential stores. Only after this third intervention, the median (effective) growth rate, λ * (t) = λ 3 − µ = −0.03 < 0 (CI [−0.05, −0.02])], finally became negative, indicating a decrease in the number of new infections. We can thus clearly relate the change points to the governmental interventions and quantify their mitigation effect. DISCUSSION We presented a Bayesian approach for a timely monitoring of the effect of governmental interventions on the spread of an epidemic outbreak. At the example of the COVID-19 outbreak in Germany, we applied this approach to infer the central epidemiological parameters and three change points in the spreading rate from the number of reported cases. We showed that change points in the spreading rate affect the confirmed case numbers with a delay of about two weeks (median reporting delay of D = 11.4 days plus a median change-point duration of 3 days). Thereby, we were able to relate the inferred change points to the three major governmental interventions in Germany: We found a clear reduction of the spreading rate related to each governmental intervention (Fig. 3), (i) the cancellation of large events with more than 1000 participants (around March 9), (ii) the closing of schools, childcare centers and the majority of stores (in effect March 16), and (iii) the contact ban and closing of all non-essential stores (in effect March 23).
Our results suggest that the full extent of governmental interventions was necessary to stop exponential growth. The first two governmental interventions brought a reduction of the growth rate λ * from 30% to 12% and down to 2%, respectively. However, these numbers still implied exponential growth. Only with the third intervention -the contact ban -we found that we have crossed the transition in new case numbers from growth to decay. However, the decay rate of about −3% (CI [−5%, −2%]) remains close to zero. Hence, even a minor increase in the spreading rate may again switch the dynamics to the unstable regime with exponential growth.
We used a formal Bayesian model comparison in order to validate the presence of change points. Our model comparison ruled out models with fewer than two change points (Table I, S2). While this may seem trivial, it has important consequences for making short-term forecasts that decision makers rely on: Demonstrating and quantifying the effect of past change points can be used to formulate priors for the effects of future, similar change points. These priors help to project the effects of more recent change points into future forecasts, even when these change points are not apparent in the reported case numbers yet. Consequently, it is important to look out for and identify potential change points as early as possible to incorporate them into forecasts.
The detection of change points and their interpretation depend crucially on an accurate estimate of the reporting delay D. Therefore, the validity of its estimate should be evaluated. In our model, D contains at least three distinct factors: the biological incubation period (median 5-6 days) [32], an additional delay from first symptoms to symptoms motivating a test (1-3 days) and a possible delay before a testing results come in (1-4 days). The sum of these delays seems compatible with our inferred median delay of D = 11.4 days, especially given the wide range of reported incubation periods.
We chose to keep our main model comparatively simple, because of the small number of data points initially available during an epidemic outbreak. With such a low number of data points, only a limited number of parameters can be effectively constrained. Hence, we chose to approximate a time-dependent spreading rate λ(t) by episodes of constant spreading rates λ i that are separated by three change points where a transition occurs. Our results show that this main model is currently sufficient for Germany: While we introduced fairly broad priors on the spreading rates, we obtained comparably narrow posterior distributions for each spreading rate λ i (Fig. 3). We additionally evaluated extensions of our main model with three change points, e.g., by explicitly taking into account the incubation period (Fig. S3). These models yield consistent results for the three change points, and all have LOO scores within one standard error of each other. Thus, we consider our main model to be sufficient for case numbers in Germany at present.
Our framework can be easily adapted to other countries and enables one to incorporate future developments. For other countries, or for forecasts within smaller communities (e.g. federal states or cities), additional details may become important, such as explicit modeling of incubation time distributions [17,37] (i.e. as done in Fig. S3), spatial heterogeneity [17,21], isolation effects [20,37], subsampling effects hiding undetected cases even beyond the reporting delay [38,39], or the age and contact structure of the population [26]. In countries where drastic changes in test coverage are expected, this will have to be included as well. The methodology presented here is capable in principle of incorporating such details. It also lends itself to modeling of continuous drifts in the spreading rate, e.g. reflecting reactions of the public to news coverage of a catastrophic situation, or people growing tired of mitigation measures. Such further adaptations, however, can only be performed on a per-country basis by experts with an intimate knowledge of the local situation. Our code provides a solid and extensible base for this. For Germany, several developments in the near future may have to be included in the model. First, people may have transiently changed their behavior over the Easter holidays; second, we expect a series of change points, as well as continuous drifts, with governments trying to ease and calibrate mitigation measures. Third, extensions to hierarchical models will enable regional assessments, e.g. on the level of federal states.
Even after the three major governmental interventions in Germany, effective growth rates remain close to zero and warrant careful consideration of future measures. At present, estimates of effective growth rates dropped to −3% and thereby remain close to zero -the watershed between exponential growth or decay. Together with the delay of approximately two weeks between infection and case report this warrants caution in lifting restrictions for two reasons: First, lifting restrictions too much will quickly lead to renewed exponential growth; second, we would be effectively blind to this worsened situation for nearly two weeks in which it will develop uninhibited. This may result in unwanted growth in case numbers beyond the level that the health system can cope with -especially when the active cases have not gone down close to zero before lifting restrictions, thus re-initiating growth from a high base level. Therefore, it is important to consider lifting restriction only when the number of active cases are so low that a two-week increase will not pose a serious threat.
In conclusion, the presented Bayesian approach allows to detect and quantify the effect of recent governmental interventions and -combined with potential subsequent interventions -to forecast future case number scenarios. Our analysis highlights the importance of precise timing and magnitude of interventions for future case numbers. It also stresses the importance of including the reporting delay D between the date of infection and the date of the confirmed case in the model. The delay D, together with the time required to implement interventions causes a total delay between an intervention and its visibility in the case numbers of about two weeks for COVID-19 in Germany. This means that changes in our behavior today can only be detected in confirmed cases in two weeks. Combined with the current spreading rate that is still around zero, the inferred spreading and observation dynamics warrant an extremely careful planning of future measures.

MATERIALS AND METHODS
As a basis for our Bayesian inference and the forecast scenarios, we use the differential equations of the well-established SIR (Susceptible-Infected-Recovered) model. We also test the robustness of our results by means of more sophisticated models, in particular an SEIR-like model that explicitly incorporates an incubation period (Fig. S3). While the SIR model-dynamics is well understood in general, here our main challenge is to estimate model parameters specifically for the COVID-19 outbreak, and to use them for forecasting. To that end, we combined a Bayesian approach -to incorporate prior knowledge -with Markov Chain Monte Carlo (MCMC) sampling -to compute the posterior distribution of the parameters and to sample from it for forecasting. Put simply, we first estimate the parameter distribution that best describes the observed situation, and then we use many samples from this parameter distribution to evolve the model equations and thus forecast future developments.
The data used comes from the Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE) dashboard [40]. The JHU CSSE provides up-to-date data on COVID-19 infections, usually a few days ahead of official German sources. The exact version of the data and code is available at [31]. Data were incorporated until April 21. Note that after this cutoff date, additional modeling of the effects of behavioral changes over the Easter holidays becomes necessary.
Simple model: SIR model with stationary spreading rate We consider a time-discrete version of the standard SIR model. In short, we assume that the disease spreads at rate λ from the infected population compartment (I) to the susceptible compartment (S), and that the infected population compartment recovers (R) at rate µ. This well-established model for disease spreading can be described by the following set of (deterministic) ordinary differential equations (see, e.g., Refs [5,6,20]). Within a population of size N , As a remark, during the onset phase of an epidemic only a very small fraction of the population is infected (I) or recovered (R), and thus S ≈ N I such that S/N ≈ 1. Therefore, the differential equation for the infected reduces to a simple linear equation, exhibiting an exponential growth Because our data set is discrete in time (∆t =1 day), we solve the above differential equations with a discrete time step (dI/dt ≈ ∆I/∆t), such that Importantly, I t models the number of all (currently) active infected people, while I new t is the number of new infections that will eventually be reported according to standard WHO convention. Importantly, we explicitly include a reporting delay D between new infections I new t and newly reported cases (C t ) as We begin our simulations at time t = 0 with I 0 infected cases and start including real-word data of reported casesĈ t from day t > D (see below for a parameterization). In our model we do not explicitly incorporate the inflow of additional infected people by travel for two reasons. First, we implicitly model the initial surge of infections with I 0 . Second, previous work showed that travel during the outbreak has only modest effects on the dynamics, e.g., travel restrictions in China merely delayed the exponential spread if not combined with reductions of spreading [41].

Full model: SIR model with weekly reporting modulation and change points in spreading rate
Our change point detection builds on a generalization of the simple SIR model with stationary spreading rate. We now assume that the spreading rate λ i , i = 1, ..., n, may change at certain time points t i from λ i−1 to λ i , linearly over a time window of ∆t i days. Thereby, we account for policy changes to reduce λ, which were implemented in Germany step by step. Thus, the parameters t i , ∆t i , and λ i are added to the parameter set of the simple model above, and the differential equations are augmented by the time-varying λ i .
In addition, we include a weekly modulation to account for lower case reports around the weekend which subsequently accumulate during the week. To model the systematic variation of case reports during the week, we adapted the newly reported cases by a reporting fraction where f w and Φ w will also be constrained by the data.

Estimating model parameters with Bayesian MCMC
We estimate the set of model parameters θ = {λ i , t i , µ, D, σ, I 0 , f w , Φ w } using Bayesian inference with Markov-chain Monte-Carlo (MCMC). The parameter σ is the scale factor for the width of the likelihood P Ĉ t θ between observed data and model (see below). Our implementation relies on the python package PyMC3 [42] with NUTS (No-U-Turn Sampling) [43] using multiple, independent Markov chains. The structure of our approach is the following: Initialization of the Markov chains via variational inference. The posterior is approximated by Gaussian distributions ignoring correlations between parameters through automatic differentiation variational inference (ADVI) [44], which is implemented in PyMC3. From this distribution, four starting points for four chains are sampled.
Burn-in phase: Each chain performs 1000 burn-in (tuning) steps using NUTS, which are not recorded. This serves as equilibration in order to sample from an equilibrium distribution in the next step.
Sampling phase: Each chain performs 4000 steps, which are used to approximate the posterior distribution. To ensure that the Chains are equilibrated and sampled from the whole posterior distribution (ergodicity), we verified that the R-hat statistic is below 1.05, which is implemented in PyMC3. The rank normalized R-hat diagnostic tests for lack of convergence by comparing the variances within chains and between chains: For identical within-chain and between-chain variances R-hat becomes 1, indicating convergence. For well-converged chains the resulting samples will describe the real-world data well, so that consistent forecasts are possible in the forecast phase.
Forecast using Monte Carlo samples. For the forecast, we take all samples from the MCMC step and continue time integration according to different forecast scenarios explained below. Note that the overall procedure yields an ensemble of forecasts -as opposed to a single forecast that would be solely based on one set of (previously optimized) parameters.
MCMC sampling details Each MCMC step requires to propose a new set of parameters θ, to generate a (fully deterministic) time series of new infected cases C(θ) = {C t (θ)} of the same length as the observed real-world datâ C = Ĉ t , and to accept or reject θ. In our case, the NUTS implementation (in PyMC3) first proposes a new set of parameters θ based on an advanced gradient-based algorithm and subsequently accepts or rejects it such that the resulting samples reflect the posterior distribution where p(Ĉ|θ) is the likelihood for the data given the parameters and p(θ) is the prior distribution of the parameters (see below). The likelihood quantifies the similarity between model outcome and the available real-world time series. Here, the likelihood is the product over local likelihoods quantifying the similarity between the model outcome for one time point t, C t (θ), and the corresponding real-world data pointĈ t . We chose the Student's t-distribution because it resembles a Gaussian distribution around the mean but features heavy tails, which make the MCMC more robust with respect to outliers [45], and thus reporting noise. The case-number-dependent width is motivated by observation noise through random subsampling [38], resulting in a variance proportional to the mean. Our likelihood neglects any noise in the dynamic process, as the SIR model is deterministic, but could be in principle extended to incorporate typical demographic noise from stochastic spreading dynamics [35,46].

Priors that constrain model parameters
As short-term forecasts are time-critical at the onset of an epidemic, the available real-world data is typically not informative enough to identify all free parameters, or to empirically find their underlying distributions. We therefore chose informative priors on initial model parameters where possible and complemented them with uninformative priors otherwise. Our choices are summarized in Tab. II for the simple model, i.e. the SIR model with stationary spreading rate for the exponential onset phase, and in Tab. III for the full model with change points, and are discussed in the following. Priors for the simple model (Table II): In order to constrain our simple model, an SIR model with stationary spreading rate for the exponential onset phase, we chose the following informative priors. Because of the ambiguity between the spreading and recovery rate in the exponential onset phase (see description of simple model), we chose a narrow log-normal prior for the recovery rate µ ∼ LogNormal(log(1/8), 0.2) with median recovery time of 8 days [20]. Note that our implementation of µ accounts for the recovery of infected people and isolation measures because it describes the duration during which a person can infect others. For the spreading rate, we assume a broad log-normal prior distribution λ ∼ LogNormal(log(0.4), 0.5) with median 0.4. This way, the prior for λ − µ has median 0.275 and the prior for the base reproduction number (R 0 = λ/µ) has median 3.2, consistent with the broad range of previous estimates [18,33,34]. In addition, we chose a log-normal prior for the reporting delay D ∼ LogNormal(log(8), 0.2) to incorporate both the incubation time between 1-14 days with median 5 [32] plus a delay from infected people waiting to contact the doctor and get tested.
The remaining model parameters are constrained by uninformative priors, in practice the Half-Cauchy distribution [47]. The half-Cauchy distribution HalfCauchy(x, β) = 2/πβ[1 + (x/β) 2 ] is essentially a flat prior from zero to O(β) with heavy tails beyond. Thereby, β merely sets the order of magnitude that should not be exceeded for a given parameter. We chose for the number of initially infected people in the model (16 days before first data point) I 0 ∼ HalfCauchy(100) assuming an order of magnitude O(100) and below. In addition, we chose the scale factor of the width of the likelihood function as σ ∼ HalfCauchy(10); this choice means that the variance in reported numbers may be up to a factor of 100 larger than the actual reported number. Priors for the full model (Table III): In order to constrain our full model, an SIR model with weekly reporting modulation and change points in the spreading rate, we chose the same priors as for the simple model but added the required priors associated with the change points. In general, we assume that each set of governmental interventions (together with the increasing awareness) leads to a reduction (and not an increase) of λ at a change point. As we cannot know yet the precise reduction factor, we adhere to assume a reduction by ≈ 50%, but always with a fairly wide uncertainty, so that in principle even an increase at the change point would be possible. We model the time dependence of λ as change points, and not as continuous changes because the policy changes were implemented in three discrete steps, which were presumably followed by the public in a timely fashion.
For the spreading rates, we chose log-normal distributed priors as in the simple model. In particular, we chose for the initial spreading rate the same prior as in the simple model, λ 0 ∼ LogNormal(log(0.4), 0.5); after the first change point λ 1 ∼ LogNormal(log(0.2), 0.5), assuming the first intervention to reduce the spreading rate by 50% from our initial estimates (λ 0 ≈ 0.4) with a broad prior distribution; after the second change point λ 2 ∼ LogNormal(log(1/8), 0.5), assuming the second intervention to reduce the spreading rate to the level of the recovery rate, which would lead to a stationary number of new infections. This corresponds approximately to a reduction of λ at the change point by 50%; and after the third change point λ 3 ∼ LogNormal(log(1/16), 0.5), assuming the third intervention to reduce the spreading rate again by 50%. With that intervention, λ 3 is smaller than the recovery rate µ, causing a decrease in new case numbers and a saturation of the cumulative number of infections.
For the timing of change points, we chose normally distributed priors. In particular, we chose t 1 ∼ Normal(2020/03/09, 3) for the first change point because on the weekend of March 8, large public events, like soccer matches or fairs, were cancelled. For the second change point, we chose t 2 ∼ Normal(2020/03/16, 1), because on March 15, the closing of schools and other educational institutions along with the closing of non-essential stores were announced and implemented on the following day. Restaurants were allowed to stay open until 6 pm. For the third change point, we chose t 3 ∼ Normal(2020/03/23, 1), because on March 23, a far-reaching contact ban (Kontaktsperre), which includes the prohibition of even small public gatherings as well as complete closing of restaurants and non-essential stores was imposed by the government authorities. Further policy changes may occur in future; however, for now, we do not include more change points.
The change points take effect over a certain time period ∆t i for which we choose ∆t i ∼ LogNormal(log(3), 0.3) with a median of 3 days over which the spreading rate changes continuously as interventions have to become effective. The precise duration of the transition has hardly any affect on the cumulative number of cases (Fig. 2 E-F). We assumed a duration of three days, because some policies were not announced at the same day for all states within Germany; moreover, the smooth transition also can absorb continuous changes in behavior.
The number of tests that are performed and reported vary regularly over the course of a week and are especially low during weekends. To account for this periodic variation, we modulated the number of inferred cases by the absolute value of a sine function with in total a period of 7 days. We chose this function as it is a non-symmetric oscillation, fitting the weekly variation of cases on a phenomenological level. For the amplitude of the modulation we chose a weakly informative Beta prior: f w ∼ Beta(mean = 0.7, std = 0.17) and for the phase a nearly flat circular distribution: Φ w ∼ vonMises(mean = 0, κ = 0.01).
Since change point detection entails evaluating models with different numbers of parameters, some form of fair model comparison is needed. This is necessary to compensate for the higher flexibility of more complex models, as this flexibility carries the risk of overfitting and overconfident forecasts. The standard approach to avoid over-fitting in machine learning is cross-validation, and cross validation has recently also been advocated for Bayesian model comparison (e.g. [3,36]), especially for models employed for predictions and forecasts. Thus, one would ideally like to compare the models with different numbers of change points by the probability they assign to previously unobserved data points. Technically this is measured by their out-of-sample prediction accuracy, i.e. their log pointwise predictive density (lppd): where the vector [y os 1 , . . . , y os 1 ] is a an out-of-sample dataset of N new data points, and where p post (θ) = p post (θ|y, M j ) is the posterior distribution of the parameters, given the in-sample data y and the model M j . In practice, the integral is approximated using a sufficient amount of samples from p post (θ). However, this approach is only reasonable if a sufficient amount of out-of-sample data is available, which is not the case in the early stages of a disease outbreak. Therefore, the pointwise out-of-sample prediction accuracy was approximated using Leave-one-out cross-validation (LOO) in PyMC3 to compute equation 6 individually for each left out data point based on the model fit to the other data points. The sum of these values, multiplied by a factor of −2 then yields the leave-one-out cross-validation (LOO-CV) score. Thus, lower LOO-CV scores imply better models. interests. Data and materials availability: Both data [48] and analysis code [31] are available online. This work is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/. This license does not apply to figures/photos/artwork or other content included in the article that is credited to a third party; obtain authorization from the rights holder before using such material.       Fig. 3 (three change points, main text) with the same model "SIR main" but with two change points. With two change points, the onset time of the second change point is close to the (middle) one inferred when using the 3-c.p. model. However, the effective growth rate λ * after the respective last c.p. is the same in both cases, λ * = −0.03 (see also Table S2) Fig. 3 (three change points, main text) but with a more involved SEIR-like model and three change points, code available online [31]. Arguably, the SEIR-like model provides are more realistic (but also more complex) description of virus propagation. It yields a slightly better (lower) LOO-score than our "SIR main" model in the cross-validation, Table S2. However, inferred parameters are compatible with the simpler model and the inferred dynamics are similar. Model details: The SEIR-like model builds on our "SIR main" model (which includes a weekend correction, see Methods). The SEIR-like model features an explicit log-normal incubation period and a lognormal reporting delay. The incubation period is implemented as a discrete convolution of multiple exposed pools with a lognormal kernel. The discrete lognormal kernel is parameterized as follows (to match the characteristic incubation time of COVID-19 [32]): median Dinc, scale parameter 0.418 and normalized to 1. The median is a free parameter with prior Normal(5, 1) (days) [32]. The reporting delay is implemented in a similar manner: as a convolution of the number of new cases with a lognormal kernel and the scale parameter is fixed to 0.3. In order to match the total delay of the main model (between the infection and the observation), the median D is a free parameter with prior LogNormal(5, 0.2) (days). Note that because of the lognormal-distributed incubation period in the SEIR-like model, the spreading and recovery rates are not directly comparable to the SIR models. We correspondingly adapted the respective priors [17]: We adapted the prior median of the recovery rate µ to 1/3 with a scale factor of 0.3, and the prior median for λ0, λ1, λ2, and λ3 to 2.0, 1.0, 0.5 and 0.25, respectively, with a scale parameter of 1 each. A: Time-dependent model estimate of the effective growth rate λ * (t). B: Comparison of daily new reported cases and the model (inset: log-lin scale). C: As B but for total (cumulative) reported cases. D-G: Prior and posterior distributions of all free parameters.  Fig. 3 (three change points, main text) with an SIR model that excludes the weekend modulation but features the three change points. In this version of the model, we excluded the assumption that daily new reported cases depend on the weekday (which is modeled as an absolute sine with an amplitude and a phase shift as inferred parameters in the main model). While the inferred parameters from the model that excludes the weekend modulation (in particular rates and onset times of change points) match the model that includes the modulation, the LOO-scores of the cross-validation are worse,   Fig. 3 (main text, three change points, same model) but with a prior for the reporting delay that is 4 times wider (panel F, third column). All parameters and change points are constrained by data. In particular, the posterior distribution of the reporting delay is slightly wider than with the original prior, but it is well constrained by data. A: Time-dependent model estimate of the effective growth rate λ * (t). B: Comparison of daily new reported cases and the model (inset: log-lin scale). C: As B but for total (cumulative) reported cases. D-G: Prior and posterior distributions of all free parameters.  Fig. 3 (main text, three change points, same model) but with a prior for the change times that is 14 days wide, instead of ∼ 2 days, (panel G, second column). All parameters and change points are constrained by data but the change points occur at later times compared to the original priors. A: Time-dependent model estimate of the effective growth rate λ * (t). B: Comparison of daily new reported cases and the model (inset: log-lin scale). C: As B but for total (cumulative) reported cases. D-G: Prior and posterior distributions of all free parameters.

Prior Posterior
FIG. S7. Sensitivity analysis: Change-point detection as in Fig. 3 (main text, three change points, same model) but with a prior for the change duration that is 4 times wider (panel G, third column). The duration of the first change point ∆t1 is robust to the wider prior; the data constrains the posterior. The durations of the second and third change point, ∆t2 and ∆t3 are not constrained by the data but they depend on our chosen priors. This is also visible in the lack of plateaus in the effective growth, panel A. However, the inferred spreading rate λ and the forecast are not sensitive to the wider priors. A: Time-dependent model estimate of the effective growth rate λ * (t). B: Comparison of daily new reported cases and the model (inset: log-lin scale). C: As B but for total (cumulative) reported cases. D-G: Prior and posterior distributions of all free parameters.