key: cord-0573189-r23zn71h
authors: Budanur, Nazmi Burak; Hof, Bjorn
title: An autonomous compartmental model for accelerating epidemics
date: 2021-10-12
journal: nan
DOI: nan
sha: d1a2666fa23187b6e983b423aa97a6e339c02cf5
doc_id: 573189
cord_uid: r23zn71h

In Fall 2020, several European countries reported rapid increases in COVID-19 cases along with growing estimates of the effective reproduction rates. Such an acceleration in epidemic spread is usually attributed to time-dependent effects, e.g. human travel, seasonal behavioral changes, mutations of the pathogen etc. In this case however the acceleration occurred when counter measures such as testing and contact tracing exceeded their capacity limit. Considering Austria as an example, here we show that this dynamics can be captured by a time-independent, i.e. autonomous, compartmental model that incorporates these capacity limits. In this model, the epidemic acceleration coincides with the exhaustion of mitigation efforts, resulting in an increasing fraction of undetected cases that drive the effective reproduction rate progressively higher. We demonstrate that standard models which does not include this effect necessarily result in a systematic underestimation of the effective reproduction rate.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1] and its variants [2] continue to challenge our lives [3] even after the rollout of several effective vaccines [4] . As of Fall 2021, non-pharmaceutical interventions such as mask mandates and travel restrictions at varying degrees still remain in place around the globe [3] . These interventions are guided by almost-real-time assessment of the ongoing epidemiological situation which rely on surveillance data and mathematical models and, are thus, prone to their uncertainties and shortcomings [5] . It is, therefore, crucial for decision making that the epidemiological models are sufficiently simple to be used in a fast-changing environment while containing the necessary amount of complexity to capture all essential features of the real epidemic.

Simple epidemic models divide a population into "compartments" according to individuals' epidemiological status and specify the rules by which the disease progresses within an individual and spreads over the population [6] . In the most basic form, these rules are given as transition rates between the compartments which can be translated into a set of ordinary differential equations (ODEs). One such model is the SEIR model where the compartments correspond to those who are susceptible (S) to infection, exposed (E) to the pathogen (but not yet a spreader), infectious (I), and removed (R) from epidemic dynamics (dead or immune). The SEIR model can capture the initial exponential increase of infections when the majority of a population is susceptible and the subsequent slow down of spreading due to the continuously increasing removed population.

In Fall 2020, the second wave of the COVID-19 in Austria exhibited a remarkably different behavior than the one that we described above. What appeared to be slowbut-steady initial increase was followed by an acceleration of the epidemic indicated by a faster-than-exponential growth in case numbers [7] . For simple mathematical models, such an observation indicates a change in rules, i.e. increased transmissibility, which can be due to seasonality [8, 9] or occurrence of more infectious variants [10] . In contrast to such expectations, here we show that the accelerating increase of COVID-19 cases reported in Austria during Fall 2020 can be captured in an autonomous compartmental model described by a time-independent (deterministic) set of ODEs. As we explain in the following, the key modeling ingredient for accelerating epidemics is explicit inclusion of mitigation efforts and their capacity limitations within the model dynamics.

Since the early stages of the pandemic, case and contact isolation has been one of the primary public health responses [11] . Using stochastic agent-based epidemic models on networks, Scarselli et al. [12] showed that while testing and contact isolation slow down an epidemic, limiting the number of available tests fundamentally alter its nature by changing the epidemic transition for large populations (in the thermodynamic limit) from gradual (second order) to sudden (first order). The driving mechanism of this qualitative change was the acceleration of transmissions when the models' testing capacity were exhausted. In the present work, we introduce the SEIRTC model which, in addition to the S, E, I, and R, includes separate compartments for the tested individuals (T) and confirmed cases (C) similar to [12] . Differently from the stochastic network models used in [12] , the dynamics of the SEIRTC model is determined by a set of ODEs, rendering it suitable for working with real data. We show that the SEIRTC model can be fitted to the COVID-19 surveillance data published by the Austrian Agency for Health and Food Safety (AGES) [13] and capture the epidemic acceleration observed in Fall 2020 without the need for a temporal modification of the infectiousness. Our results suggests that during this period, the effective reproduction rate, i.e. the average number of secondary cases originating from a primary one were systematically underestimated.

Before introducing the SEIRTC model, we briefly recapitulate the standard SEIR model represented by the state transition diagram in Figure 1A . The corresponding ODEs can be written by expressing the rates of changes in compartments' populations according Figure 1 : State transition diagrams of SEIR (A) and SEIRTC (B) models where the encircled letters denote the compartments into which the population is divided and the arrows along with their labels underneath indicate the transition rates between the compartments. The compartments are S: susceptible, E: exposed, I: infectious, R: removed (recovered or dead), I a/s : a/symptomatic infectious, R u/k : un/known removed, T S/E/I : tested susceptible/exposed/infectious, C: case.

to the transition rates implied by the annotated arrows in Figure 1A aṡ

where,˙denotes derivative with respect to time, N = S + E + I + R is the population, β is the transmission parameter, and γ E and γ I are the inverse latent and infectious times, respectively. The transmission parameter β can be interpreted as the number of interactions per person per unit time multiplied by the transmission probability at each interaction [6] . The underlying assumption of SEIR model is that the population is well mixed, thus, the effects solely due to the network heterogeneities are neglected. The SEIR model forms the basis of our model depicted in Figure 1B where we additionally have test (T) and case (C) compartments. Moreover, I, T , and R compartments are split into sub-compartments as follows. I s and I a refer to the "symptomatic" and "asymptomatic" infectious individuals, respectively, the latter of whom are those who show no symptoms throughout their infectious period and spread the disease at a relative risk ρ as implied by the S → E term. The subscripts of the test sub-compartments indicate the epidemiological state of the tested individuals, thus, the susceptibles who are tested (T S ) return to S, whereas the exposed and infectious individuals become "cases" (C) after being tested. Finally, the removed individuals are split into "known" (R k ) and "unknown" (R u ) parts for convenience in presenting our results to follow.

Similar to those in the SEIR model, the parameters γ x of the SEIRTC model refer to the inverse mean lifetime at the compartment x, with the exception γ s which is the inverse mean symptom onset time, at which point a symptomatic infectious individual becomes detectable. Hence, this delay term accounts for presymptomatic transmissions which are believed to play a significant role in spreading COVID-19 [14] . γ T is the test turn-over time, i.e. the time from the administration of a test to its outcome. Finally, γ C is not an independent parameter but is determined as

where we assume that at time t, the ratio of cases that are identified before and after developing symptoms are proportional to the number of individuals in T E and T I , respectively. In doing so, we neglect a delay by γ −1 T and avoid working with delaydifferential equations.

The transition rates in Figure 1B have several probabilistic factors. These are p: probability of an asymptomatic infection, ρ: relative risk of transmission from an asymptomatic individual, d: probability of testing a symptomatic infectious individual, g: probability of testing an exposed individual before becoming infectious. All but the last one of these probabilities are independent parameters to be determined via literature estimates or model fitting. Because the detection of an infection before developing any symptoms can only be possible via contact tracing, the probability g is a complex function of the number of (un)identified infections, the underlying social network structure, and the contact tracing capacity. We propose the ansatz

where κ and T m are fit parameters. The fraction C/(C +I s +I a ), i.e. the ratio of known cases to the total number of infectious individuals, appears in (3) multiplicatively since the larger this ratio is the higher becomes the probability that a new case is originated from a known one. Finally, the term in square brackets models the contact tracing capacity as it is 1 for T m T S + T E and 0 for T m T S + T E . Similar to g, the number of susceptible individuals to be tested per unit time f is also an unknown function of the social network and contact tracing procedures. As our second ansatz, we take

where the term αS(C/N ) is analogous to the S → E term of the SEIR model (Figure 1A) , thus, we interpret the parameter α as the number of interactions per person per day multiplied by the probability of no infection to occur and the probability of contact identification when T m T S + T E . The limiting term in the square brackets is identical to that in (3) .

With the ansätze (3) and (4), the SEIRTC ODEs similar to (1) corresponding to the state transition diagram Figure 1B can be obtained by expressing the rates of change in compartment populations as the shown transition rates. In our numerical results to follow, we simulate these using odeint function of scipy [15] .

In a real-life scenario, testing constitutes the primary source of information as most countries publish the daily numbers of tests they conduct and those with a positive outcome (incidence). For model fitting, we make use of both of these measurements and determine the initial state of the population and free parameters through a weighted least squares fit [16] that minimizes the cost function

where [n] denotes the discrete time in days,˜indicates measurements coming from the surveillance data, w T,C are weights, and C

is the number of new cases recorded in the model on day n. In order to reduce the number of fit parameters, we make the following simplifying assumptions. Whenever available, we take literature values for parameters or restrict them to the established estimated interval. While the number of exposed individuals on day 0 is varied, the initial populations of T S , T E , I a , I s , T I , C are adjusted through a fixed-point iteration that transfers individuals from S to these compartments while minimizing the error between the simulated dynamics of the first 10 days from exponential fits. R k (0) and R u (0) are both set to the number of cases registered until the first day. Although this is an arbitrary assumption for R u , it has no effect on the dynamics as long as it is much less than the total population, which is the case for our results to follow. With these assumptions, the set of fit parameters becomes {E(0), α, κ, T m , β, γ s , d}.

In order to estimate the effective reproduction number (R t ) from surveillance and model data we utilize the python implementation [17] of the Cori et al. [18] 's EpiEstim algorithm that is based on the Bayesian inference of R t from a Gamma-distributed prior assuming Poisson-distributed transmissions. This method was also used by AGES [19] who performs the real-time epidemilogical monitoring of the ongoing COVID-19 situation in Austria.

We consider the second wave of COVID-19 in Austria from September 1 to November 3, 2020, on which day the country went into its second lockdown in order to protect its healthcare system from an otherwise-inevitable overload. Figure 2A , B shows the 7-day moving averages of the numbers of confirmed cases and performed tests during this period, respectively (retrived from [13] ). Fits to these data by the SEIRTC model with parameters in Table 1 are also shown in Figure 2A Figure 2A , we also show fits by SEIR models with N , γ E , and γ I as in Table 1 and w T [n] = 0, w C [n] = 1 (SEIR 2 ). When unit weights are chosen, the SEIR model with the best fit parameters E 2 (0) = 535 , β 2 = 0.332, underestimates the initial case numbers whereas when the weights are inversely proportional to the daily number of confirmed cases, the later case numbers are underestimated with the best fit parameters E 1 (0) = 1306 , β 1 = 0.282. In contrast, the fit by the SEIRTC model captures both episodes. As shown in Figure 2B , the SEIRTC model also captures the reported test numbers with a visible change in its trend during the days following September 15. This coincides with the exhaustion of the contact tracing capacity as reflected by the plateaus of T S and T E in Figure 2C . Consequently, the ratio I [n]/C [n] (Figure 2D) 

, of the daily new infections to those that are detected starts growing, yielding more and more infections that remain undetected. An observable consequence of this is the increase of the proportion of the tests with a positive outcome as shown in Figure 2E where we plotted this quantity using the data reported and those obtained from the SEIRTC model.

In Figure 2F , we show the estimates of the effective reproduction number R t calculated using the case data shown in Figure 2A and those computed from the SEIRTC model's daily case (orange, dashed) and infections (green, dotted-dashed). In these calculations, we chose a smoothing window of 7 days and used a serial interval obtained by discretizing a Gamma distribution with a mean 4.46 and a standard deviation 2.63 days as estimated by AGES [19] . As a simple model of a nonpharmaceutical intervention, we investigated hypothetical lockdown scenarios during which the transmission parameter β is reduced by half. Specifically, we started SEIR 2 and SEIRTC models with their initial conditions and model parameters same as those that we used in Figure 2 , initiated lockdowns when the 7-day incidence (per 100,000 people) reached a certain value, and measured the lockdown duration necessary for reducing the incidence by 10. The results of these simulations are shown in Figure 3 . Figure 2A , B demonstrates that the SEIRTC model with parameters listed in Table 1 and the initial population shown in Figure 2C is able to capture the accelerating spread of COVID-19 observed in Austria in Fall 2020. We note that some of these fit parameters are correlated, hence, it is possible to find different set of parameters that result in comparable fitting errors (5) . For example, one can easily imagine increasing Figure 3 : Number of days in "lockdown" necessary for reducing the 7-day incidence (per 100,000) by 10 as a function of the incidence at which the lockdown begins. Lockdown is modeled by a drop of the transmission parameter β to half of its value.

the ratio d of the symptomatic cases that are being detected while reducing the time constant γ s for symptomatic case detection and having a similar-quality fit to the reported data. For predicting the future of an ongoing outbreak, quantification of such parameter correlations and uncertainties is crucial and, thus, should be carefully taken into consideration [24] . Because near-real-time prediction is not our goal here, we focus on the qualitative aspects of our findings that are robust to parameter uncertainties.

While Figure 2A illustrates how an autonomous SEIR model cannot capture an accelerating epidemic, it also suggest how a nonautonomous SEIR model with a timevarying transmission parameter β(t) could have indeed describe the observed case numbers. One could then have interpreted the accelerating spread as being due to the seasonal effects such as people spending more time indoors hence increasing chance of transmission. Another -somewhat trivial -modification to the SEIR model could have been addition of infectious individuals to the model by hand as a proxy for people bringing the virus from outside the country through travel as suggested by the Austrian then-Chancellor Sebastian Kurz who claimed that the sudden increase of the COVID-19 cases during Fall 2020 was largely due to the Austrians of foreign origin who brought the virus back from their their countries of origin [25] . It is conceivable that both of these factors have played some role in the sudden spread of COVID-19 in Austria in Fall 2020 however they as such do not explain the coinciding increase in the test-positive rate. As we point out in this study capacity limits in mitigation have played a key role in the epidemic acceleration and offer an explanation for the observed events including the change in positive rate.

The methods for estimating the effective reproduction number are usually believed to be robust against incomplete observations under the assumption that an approximately same fraction of infections are recorded on consecutive days [26] . This is also reflected in our initial estimates of the effective reproduction number, since when the ratio of the number of infections to that of confirmed cases is constant as in the ini-tial phase of Figure 2D , R (C) t and R (I+C) t coincide in Figure 2F 1 . This, however, breaks down as soon as the contract tracing limit is reached; after this point, the effective reproduction number based only on the case data systematically underestimates the one that is based on the actual number of infections. Although the difference between R (C) t and R (I+C) t appear small in Figure 2F , it should be recalled that this difference translates into the number of infections exponentially, thus, has a dramatic real-life consequences such as expected number of hospital admissions.

As one should expect, the increase of undetected cases in Figure 2D coincide with that of the ratio of positive tests in Figure 2E , which is observable during an epidemic. While this information could, and probably should, be incorporated into statistical methods for estimating R t , we believe that an increasing positive test ratio is a sufficient reason for a dramatic intervention such as a lockdown, which is essentially inevitable once the contact-tracing capacity is exhausted.

During the second wave of covid-19 in Austria, the policy makers insisted that a lockdown would be the last option in the country's pandemic response [27] . Decreasing lockdown durations for the SEIR model as shown in Figure 3 might indeed suggest this as a reasonable compromise to minimize the number of days during which the economic and social activities are halted. Note, however, that this behavior changes dramatically in the SEIRTC model since the uncontrolled spread following the breakdown of contact tracing makes it progressively harder to reduce the case numbers. This is why we believe that a steadily increasing ratio of positive tests necessitates a lockdown. At that point, early action not only saves lives but also shortens the lockdown duration necessary to regain control.

SARS-CoV-2 (COVID-19) by the numbers. eLife, 9

World Health Organization. Tracking SARS-CoV-2 variants. URL https:// covid19-dashboard.ages.at. Retrieved on 25 Augusut

A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker)

A global database of COVID-19 vaccinations

The Uncertain Future in How a Virus Spreads

Modeling Infectious Diseases in Humans and Animals

Woran es liegt, dass die corona-zahlen so schnell anstiegen

Seasonality and uncertainty in global COVID-19 growth rates

The role of seasonality in the spread of COVID-19 pandemic

Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G increases Infectivity of the COVID-19 Virus

Feasibility of controlling CovId-19 outbreaks by isolation of cases and contacts

Discontinuous epidemic transition due to limited testing

AGES Dashboard COVID-19

Presymptomatic SARS-CoV-2 infections and transmission in a skilled nursing facility

SciPy 1.0: fundamental algorithms for scientific computing in python

Modeling and Inverse Problems in the Presence of Uncertainty

Estimation and worldwide monitoring of the effective reproductive number of SARS-CoV-2

A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics

Epidemiologische Parameter des COVID-19 Ausbruchs -Update

Abgestimmte Erwerbsstatistik und Arbeitsstättenzählung

Asymptomatic SARS-CoV-2 infections: a living systematic review and meta-analysis

Modeling shield immunity to reduce COVID-19 epidemic spread

Fragen zur COVID-19 testung -stand

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts

Kurz blames austrians 'of foreign origin' for high infection rates

Practical considerations for measuring the effective reproductive number, R t

Chronology of the corona crisis in austria -part 4: Lockdowns, mass testing and the launch of the vaccination campaign