Delays, Masks, the Elderly, and Schools: First Covid-19 Wave in the Czech Republic

Abstract

Running across the globe for nearly 2 years, the Covid-19 pandemic keeps demonstrating its strength. Despite a lot of understanding, uncertainty regarding the efficiency of interventions still persists. We developed an age-structured epidemic model parameterized with epidemiological and sociological data for the first Covid-19 wave in the Czech Republic and found that (1) starting the spring 2020 lockdown 4 days earlier might prevent half of the confirmed cases by the end of lockdown period, (2) personal protective measures such as face masks appear more effective than just a realized reduction in social contacts, (3) the strategy of sheltering just the elderly is not at all effective, and (4) leaving schools open is a risky strategy. Despite vaccination programs, evidence-based choice and timing of non-pharmaceutical interventions remains an effective weapon against the Covid-19 pandemic.

1 Introduction

The Covid-19 pandemic of the SARS-CoV-2 virus has held the world in its grip for nearly 2 years. With the first three cases reported on March 1, 2020, the epidemic in the Czech Republic (Czechia) started and was initially fueled by Czech citizens returning from the alpine ski resorts of Italy and Austria. Population-wide interventions began on March 11, 2020, with the closing of schools, soon followed by travel restrictions, closing of restaurants, sports and cultural facilities and shops (with some exceptions), as well as the introduction of a duty to wear face masks and keep at least 2 m inter-personal distance in public (Table 1). By May 2020 the epidemic situation in Czechia stabilized and the lockdown restrictions were gradually relaxed. During the summer months only sporadic local outbreaks occurred and life returned nearly to normal. Since then, other four waves vexed Czechia, resulting in unprecedented numbers of hospital admissions and deaths (Ministry of Health of the Czech Republic 2020; Komenda et al. 2020).

Table 1 Population-wide interventions against Covid-19 (lockdown) in the Czech Republic during the first (spring 2020) wave

Since lockdowns not only damage national economy but also negatively affect the mental state of people (Pierce et al. 2020), governments all over the world tend not to implement severe population-wide restrictions until they become unavoidable. Fortunately, we can learn from past epidemic waves and increase the effectiveness of any forthcoming interventions. How fast do we need to react in order to limit potential cases and casualties? How can we protect ourselves effectively or limit social contacts when wearing hazmat suits or locking ourselves at home for 2 or 3 weeks are not viable options? Can we limit the worst impacts of the epidemic just by sheltering the elderly with the rest living more or less normally? What about leaving schools open or not forcing people to work from home? These questions have long sprouted controversies among politicians and in the general public not only in Czechia, and appear to have not been explored deeply in the scientific literature.

Here, we develop a mathematical model of the Covid-19 epidemic to address these four questions. The model is structured by age and type of inter-individual contacts (at home, school, work and in the community), and considers all important epidemic classes. We parameterize the model by combining the public health data on the first (spring 2020) wave in Czechia, including the type and timing of interventions implemented, sociological data from surveys carried out before and during this wave, and data published in the scientific literature. The major reason to use this wave is a relatively clear start and stability of both epidemic and behavioral characteristics: population-wide interventions were implemented soon after its outbreak and stayed virtually unchanged until its decline, a situation that has not echoed since in such clarity.

In order to account for all these sources of information and their inevitable uncertainty we used the Approximate Bayesian Computation (ABC) framework. This technique, used to estimate parameters of complex models in genomics and other biological disciplines (Beaumont 2010; Csilléry et al. 2010), including epidemiology (Toni et al. 2009; Blum and Tran 2010; Luciani et al. 2009), allows to assess uncertainty in model simulations.

2 Methods

Our epidemic model is an extension of the classic SEIR model, structured by age and type of contacts. We start with describing an unstructured version of the model. Its extension to age and various types of contacts then follows. A general scheme of our model is provided in Fig. 1.

Fig. 1

A general scheme of the epidemic model developed in this study

Core model Due to contacts with infectious individuals, susceptible individuals (S) may become exposed (E), that is, infected but not yet infectious (the process of infection transmission is described below). The exposed individuals then become either asymptomatic for the whole course of infection ($I_{\mathrm{a}}$, with a probability $1-p_{\mathrm{S}}$) or presymptomatic for just a short period of time before becoming symptomatic ($I_{\mathrm{p}}$, with a probability $p_{\mathrm{S}}$). The $I_{\mathrm{p}}$ individuals later become symptomatic. We assume that a proportion $p_{\mathrm{T}}$ of such symptomatic individuals decide to undergo testing for the presence of SARS-CoV-2 ($I_{\mathrm{s}}$). On the contrary, the proportion $1-p_{\mathrm{T}}$ of symptomatic individuals (most likely those with very mild symptoms) decide not to undergo testing and stay at home ($I_{\mathrm{h}}$). Since in the Czech Republic, deaths attributed to Covid-19 did not generally occur outside hospitals, the $I_{\mathrm{a}}$ and $I_{\mathrm{h}}$ individuals eventually recover (R).

Considering discrete time, with one time step corresponding to 1 day, our model so far consists of the following equations:

$$\begin{aligned} \begin{array}{l} S[t+1] = S[t] - \lambda \,S[t] - L[t-8]/p_{\mathrm{S}}, \\ E[t+1] = E[t] + \lambda \,S[t] - \sigma \,E[t] + L[t-8]/p_{\mathrm{S}}, \\ I_{\mathrm{a}}[t+1] = I_{\mathrm{a}}[t] + (1-p_{\mathrm{S}})\,\sigma \,E[t] - \gamma _{\mathrm{a}}\,I_{\mathrm{a}}[t], \\ I_{\mathrm{p}}[t+1] = I_{\mathrm{p}}[t] + p_{\mathrm{S}}\,\sigma \,E[t] - \xi \,I_{\mathrm{p}}[t], \\ I_{\mathrm{h}}[t+1] = I_{\mathrm{h}}[t] + (1-p_{\mathrm{T}})\,\xi \,I_{\mathrm{p}}[t] - \gamma _{\mathrm{s}}\,I_{\mathrm{h}}[t], \\ I_{\mathrm{s}}[t+1] = I_{\mathrm{s}}[t] + p_{\mathrm{T}}\,\xi \,I_{\mathrm{p}}[t], \\ R[t+1] = R[t] + \gamma _{\mathrm{a}}\,I_{\mathrm{a}}[t] + \gamma _{\mathrm{s}}\,I_{\mathrm{h}}[t], \end{array} \end{aligned}$$

(1)

where the parameters $\sigma$, $\xi$, $\gamma _{\mathrm{a}}$, and $\gamma _{\mathrm{s}}$ represent probabilities with which individuals leave the respective model classes. All model variables are summarized in Table 2, and all model parameters in Tables 3 and 4.

Table 2 List of state and testing-related variables used in the model

The force of infection $\lambda$ in model (1) sums contributions from all infectious classes, that is, $I_{\mathrm{a}}$, $I_{\mathrm{p}}$, $I_{\mathrm{h}}$, and $I_{\mathrm{s}}$:

$$\begin{aligned} \lambda = \beta \,C\,\frac{r_{\beta }\,I_{\mathrm{a}}[t] + r_{\beta }\,I_{\mathrm{p}}[t] + r_{\mathrm{C}}\,I_{\mathrm{h}}[t] + r_{\mathrm{C}}\,I_{\mathrm{s}}[t]}{N}. \end{aligned}$$

(2)

Here, $\beta$ is a probability of infection transmission upon contact between susceptible and infectious individuals, C is the contact rate (the mean number of other individuals an individual has an effective contact with per day), $r_{\beta }$ is a factor reducing the infection transmission probability for an asymptomatic individual relative to a symptomatic one, $r_{\mathrm{C}}$ is a factor reducing the contact rate of a symptomatic individual relative to an asymptomatic one (having symptoms should force an individual to reduce contacts with others), and N is the total population size. Positively tested individuals that are isolated, at home or in hospitals (see below), and those that die, should not be counted in N. However, since numbers of such individuals are at any time are negligible relative to the population size in Czechia, we disregard this complexity (results are not at all affected by this assumption).

Initial state The hitherto unexplained variable L[t] in model (1) accounts for the imported Covid-19 cases from abroad, mostly from Italy and Austria. A list of all confirmed (symptomatic) imported cases is available at http://onemocneni-aktualne.mzcr.cz/covid-19. However, we do not introduce such imported cases as symptomatic. Rather, we assume they came earlier as exposed, and introduce them before they were actually tested positive (to account for a delay between exposition and confirmation). Moreover, to account for the likely situation that some of the imported cases remained undetected as being asymptomatic for the whole course of infection, we divide the number of known imported cased by $p_{\mathrm{S}}$, the probability of exposed individuals eventually becoming symptomatic (Table 5).

Table 3 List of model parameters

Observational layer This layer models delays in reporting infectious individuals. The period from the onset of symptoms, through sampling and subsequent processing, up to infection confirmation and case isolation is assumed to take $d_{\mathrm{T}}$ days (testing delay). We consider all reported infections as symptomatic, since during the study period asymptomatic and symptomatic individuals were not distinguished in the public health data and the latter apparently dominated the total numbers. We thus redefine equation for $I_{\mathrm{s}}$ as:

$$\begin{aligned} \begin{array}{l} I_{\mathrm{s}}[t+1] = I_{\mathrm{s}}[t] + p_{\mathrm{T}}\,\xi \,I_{\mathrm{p}}[t] - \eta [t]\,I_{\mathrm{s}}[t], \end{array} \end{aligned}$$

(3)

where $\eta = 1-\exp (-1/d_{\mathrm{T}})$ is the testing rate. The testing rate $\eta [t]$ commonly increases in time, as the whole testing process becomes more efficient as the epidemic unfolds. We use data on the Covid-19 epidemic in Czechia to quantify the testing delay separately for March and for April and May (Table 3).

The number of newly (positively) tested individuals at day t thus equals $\eta [t]\,I_{\mathrm{s}}[t]$. Therefore, the total number of such individuals yet to be reported (B) is

$$\begin{aligned} B[t+1] = B[t] + \eta [t]\,I_{\mathrm{s}}[t] - \kappa \,B[t]. \end{aligned}$$

(4)

Here $\kappa$ is the publication rate, calculated as $\kappa = 1-\exp (-1/d_{\mathrm{P}})$, where $d_{\mathrm{P}}$ is the period from case confirmation to case reporting (publication delay). Although this rate may change in time, too, it was relatively constant during the spring 2020 wave in Czechia. The total number of confirmed cases reported until and including time $t+1$ (K) is therefore

$$\begin{aligned} K[t+1] = K[t] + \kappa \,B[t]. \end{aligned}$$

(5)

Dynamics in hospitals This part of the model neither enters the calibration procedure nor affects the published results, but without it the model would clearly be incomplete. A proportion $p_{\mathrm{H}}$ of (positively) tested individuals (those with relatively severe symptoms) require hospitalization, whereas the remaining proportion $1-p_{\mathrm{H}}$ are sent home to isolation ($I_{\mathrm{z}}$) until recovery. Hospitalized individuals may follow several pathways, depending on the number of hospital states one considers. Whereas many studies did not consider any class of hospitalized individuals (Kissler et al. 2020), the others considered one to three hospital states: one to cover all hospitalized individuals (Giordano et al. 2020), two to distinguish between common hospital beds and ICUs (Domenico et al. 2020), and three to further detach ICU patients that need lung ventilators or ECMOs (Weissman et al. 2020). We consider here one hospital state, for which we introduce three parameters: the proportion $p_{\mathrm{D}}$ of hospitalized individuals that eventually die, mean duration from hospital admission to death $d_{{\mathrm{HD}}}$, and mean duration from hospital admission to recovery $d_{\mathrm{HR}}$. We thus introduce two hospital classes, $H_{\mathrm{D}}$ and $H_{\mathrm{R}}$, representing the hospitalized individuals that eventually die or recover, respectively, and extend the above model equations as follows:

$$\begin{aligned} \begin{array}{l} I_{\mathrm{z}}[t+1] = I_{\mathrm{z}}[t] + (1-p_{\mathrm{H}})\,\eta [t]\,I_{\mathrm{s}}[t] - \gamma _{\mathrm{s}}\,I_{\mathrm{z}}[t], \\ [0.5ex] H_{\mathrm{D}}[t+1] = H_{\mathrm{D}}[t] + p_{\mathrm{D}}\,p_{\mathrm{H}}\,\eta [t]\,I_{\mathrm{s}}[t] - \alpha _{\mathrm{HD}}\,H_{\mathrm{D}}[t], \\ [0.5ex] H_{\mathrm{R}}[t+1] = H_{\mathrm{R}}[t] + (1-p_{\mathrm{D}})\,p_{\mathrm{H}}\,\eta [t]\,I_{\mathrm{s}}[t] - \alpha _{\mathrm{HR}}\,H_{\mathrm{R}}[t], \\ [0.5ex] R[t+1] = R[t] + \gamma _{\mathrm{a}}\,I_{\mathrm{a}}[t] + \gamma _{\mathrm{s}}\,I_{\mathrm{h}}[t] + \gamma _{\mathrm{s}}\,I_{\mathrm{z}}[t] + \alpha _{\mathrm{HR}}\,H_{\mathrm{R}}[t], \\ [0.5ex] D[t+1] = D[t] + \alpha _{\mathrm{HD}}\,H_{\mathrm{D}}[t], \end{array} \end{aligned}$$

(6)

where $\alpha _{\mathrm{HX}} = 1-\exp (-1/d_{\mathrm{HX}})$ for $X=D,R$.

Age structure As SARS-CoV-2 is known to differentially impact children, adults and seniors (Davies et al. 2020), we distinguish three major age classes: 0–19 years (children), 20–64 years (adults), and 65+ years (seniors). These classes interact via the force of infection. Both the probability of infection transmission upon contact $\beta$ and the daily number of contacts C are now $3 \times 3$ matrices, referred to below as the transmission matrix and the contact matrix, respectively.

The transmission matrix $\beta$ is assumed to have the following structure:

$$\begin{aligned} \beta = \left[ \begin{array}{lll} \beta _1 &{}\quad \beta _2 &{}\quad \beta _2 \\ \beta _2 &{}\quad \beta _3 &{}\quad \beta _4 \\ \beta _2 &{}\quad \beta _4 &{}\quad \beta _5 \end{array} \right] , \end{aligned}$$

(7)

where $\beta _1$ is a transmission probability between two children, $\beta _2$ is a transmission probability between children and adults or seniors, $\beta _3$ is a transmission probability between two adults, $\beta _4$ is a transmission probability between seniors and adults, and $\beta _5$ is a transmission probability between seniors. We estimate parameters $\beta _1$, $\beta _2$, $\beta _3$, $\beta _4$, and $\beta _5$ by fitting our model to data on the age-specific cumulative numbers of confirmed cases.

The contact matrix C describes the mean number of other individuals of any age cohort (rows) that an individual of an age cohort (columns) meets per day; Prem et al. (2017) published such a contact matrix for 152 countries, including the Czech Republic. Moreover, they expressed it as a sum of four specific contact matrices describing daily numbers of contacts at home ($C_{\mathrm{H}}$), school ($C_{\mathrm{S}}$), work ($C_{\mathrm{W}}$), and of other types of contact ($C_{\mathrm{O}}$):

$$\begin{aligned}&C_{\mathrm{H}} = \left[ \begin{array}{lll} 1.52 &{}\quad 0.67 &{}\quad 0.036 \\ 2.84 &{}\quad 2.05 &{}\quad 0.20 \\ 0.93 &{}\quad 0.58 &{}\quad 0.75 \end{array} \right] ,\ C_{\mathrm{S}} = \left[ \begin{array}{lll} 4.77 &{}\quad 0.20 &{}\quad 0.0014 \\ 1.81 &{}\quad 0.33 &{}\quad 0.0075 \\ 0.022 &{}\quad 0.019 &{}\quad 0.022 \end{array} \right] , \end{aligned}$$

(8)

$$\begin{aligned}&C_{\mathrm{W}} = \left[ \begin{array}{lll} 0.085 &{}\quad 0.19 &{}\quad 1.4\times 10^{-5} \\ 0.42 &{}\quad 5.28 &{}\quad 9.4\times 10^{-5} \\ 1.75\times 10^{-5} &{}\quad 0.00012 &{}\quad 4\times 10^{-5} \end{array} \right] ,\ C_{\mathrm{O}} = \left[ \begin{array}{lll} 1.61 &{}\quad 0.78 &{}\quad 0.24 \\ 1.10 &{}\quad 3.94 &{}\quad 1.01 \\ 0.15 &{}\quad 0.89 &{}\quad 0.93 \end{array} \right] .\nonumber \\ \end{aligned}$$

(9)

We exploit this division when defining and exploring various realistic intervention strategies.

Once infected, individuals of each age cohort proceed independently of individuals of other age cohorts. However, some model parameters that decide on specific pathways through the model are made age-specific. These are probabilities of becoming symptomatic $p_{\mathrm{S}}$, undergoing testing $p_{\mathrm{T}}$, requiring hospitalization $p_{\mathrm{H}}$, and dying $p_{\mathrm{D}}$ (Tables 3 and 4). Last but not least, since we know age of any imported Covid-19 case in the initial phase of epidemic, we assign each such case to the appropriate age cohort.

Table 4 List of parameters estimated by the ABC model calibration technique, including their corresponding prior distributions

Table 5 List of examined scenarios and the corresponding behavior-based parameters

Sociological data Our baseline scenario considers all interventions that were in effect during the lockdown initiated in March 2020 (Table 1). In modeling those interventions, we exploit a division of the contact matrix C into four matrices describing contacts at home ($C_{\mathrm{H}}$), school ($C_{\mathrm{S}}$), work ($C_{\mathrm{W}}$), and other types of contacts ($C_{\mathrm{O}}$). Starting from the corresponding dates listed in Table 1, we multiply the respective matrices by factors $r_{\mathrm{H}}=0.44$ for home, $r_{\mathrm{W}}=0.45$ for work, $r_{\mathrm{O}}=0.35$ for other types of contact, and $r_{\mathrm{S}}=0$ for schools. Moreover, personal protection, activated on March 19, 2020, including wearing face masks on public, wide use of disinfection and keeping inter-individual distances of more than 2 m on public, was modeled as follows. Denoting by $r_{\mathrm{P}}$ compliance of using personal protection measures (averaging masks and hygiene and computing mean over the high efficiency data), the chance that two randomly selected individuals do not both protect themselves is $(1-r_{\mathrm{P}})^2$, that one is protected while the other is not is $2\, r_{\mathrm{P}}\, (1-r_{\mathrm{P}})$, and that both are protected is $r_{\mathrm{P}}^2$. Assuming that infection transmission due to personal protection is in any individual reduced by a factor $p_{\mathrm{P}}$, all elements of the transmission matrix $\beta$ are then multiplied by a factor

$$\begin{aligned} r_{\mathrm{P}}^2\, p_{\mathrm{P}}^2 + 2\, r_{\mathrm{P}}\, (1-r_{\mathrm{P}})\, p_{\mathrm{P}} + (1-r_{\mathrm{P}})^2. \end{aligned}$$

(10)

We use $r_{\mathrm{P}}=0.88$ that represents 88% compliance of using personal protection measures and estimate $p_{\mathrm{P}}$ via model calibration (Table 4). These numbers except schools are based on results of public opinion surveys organized by two agencies during lockdown: panel surveys by the PAQ Research agency (www.paqresearch.cz) and retrospective questioning by the Median agency (www.median.eu). Results of the former survey are summarized in Fig. 2. This data show that during the second half of March and essentially during the whole April, contacts of all kinds have been largely reduced while personal protection was significant.

Fig. 2

Results of a sociological survey we use in our study. Behavioral responses before and during the lockdown in Czechia, based on panel surveys organized by the PAQ Research agency (www.paqresearch.cz) every week (in between the subsequent surveys, data were linearly interpolated). The orange line is a proportional reduction in the numbers of social contacts relative to the pre-pandemic state (value 1), the red and blue lines are proportions of respondents that reported using face masks and increased hygiene (Color figure online)

Model calibration Values of several model parameters remain uncertain, of which the transmission matrix $\beta$ is virtually always one of them. This and some other model parameters, listed in Table 4, are estimated by fitting the simulated cumulative number of confirmed cases in each age class to the age-specific time series on the actually reported cumulative numbers of confirmed cases in the Czechia.

There are many ways to meaningfully perform model calibration on real-world data and many optimization and filtering methods exist (Yang et al. 2014). We use the Approximate Bayesian Computation (ABC), a technique used to estimate parameters of complex models in genomics and other biological disciplines (Toni et al. 2009; Beaumont 2010; Csilléry et al. 2010; Blum and Tran 2010; Luciani et al. 2009). The major advantage of this method is that it naturally works with all sources of uncertainty acknowledged in the model. At the same time, the ABC does not rely on likelihood calculations and in case of sufficient computation power it can be used with models of virtually any complexity. The variant of ABC with rejection sampling that we use consists of three steps. First, we use our model to simulate summary statistics for calibration (the age-specific numbers of confirmed cases; K) 200,000 times, drawing the uncertain model parameters from prior distributions based on literature and available data on the Czech Republic epidemic; selected prior distributions for the parameters to be estimated are given in Table 4. Second, we compare the simulated summary statistics with the observed one, using the Euclidean distance D. Third, we select model simulations that satisfy $D < \epsilon$, where $\epsilon$ was chosen to pass 0.05% (100) of simulations into the selected set. Given that the used summary statistics are informative, the distribution of parameters corresponding to the selected simulations is known to converge from outside to the Bayesian posterior distribution of parameter values with $\epsilon$ going to 0, and is referred to as the approximate posterior (Beaumont 2010). The choice of $\epsilon$ and the number of simulations in the ABC is driven by compromise between computation power and smoothness and accuracy of the approximate posterior.

The set of selected parametric sets thus allows us to evaluate remaining parameter uncertainty, given the available data and adopted summary statistics (Toni et al. 2009; Beaumont 2010; Csilléry et al. 2010). This is crucial to realize, since although different parameter sets may similarly fit the available data (have similar summary statistics), and often provide similar short-term predictions, they may demonstrate significant differences in longer term and in interplay with intervention policies. To cope with such uncertainty, we do not evaluate only an absolute impact of an intervention scenario, but calculate also its impact relative to the baseline scenario, separately for each selected parametric set. If that relative impact is consistent over the whole posterior distribution of parameters, we may have confidence in its potential effect on the epidemic spread. To apply the ABC technique, we use the abc package in R (Csilléry et al. 2015), modified to work with non-normalized summary statistics.

Parameterizing individual scenarios In scenarios R1 (all measures set 4 days earlier than in reality) and R2 (all measures set 4 days later than in reality), the behavior-based parameters $r_{\mathrm{H}}$, $r_{\mathrm{S}}$, $r_{\mathrm{W}}$, $r_{\mathrm{O}}$, and $r_{\mathrm{P}}$ stayed the same as in the baseline scenario, yet all dates in Table 1 were shifted either 4 days earlier (R1) or 4 days later (R2). On the contrary, the remaining scenarios R3–R8 keep the dates of measures’ initiation as in the baseline scenario, but instead the parameters $r_{\mathrm{H}}$, $r_{\mathrm{S}}$, $r_{\mathrm{W}}$, $r_{\mathrm{O}}$, and $r_{\mathrm{P}}$ vary. All scenarios are described in detail in Table 5.

3 Results

3.1 Using Initial Phase and Lockdown Period as Baseline Scenario

Our baseline scenario, to which we compare a number of alternatives, concerns the initial phase of the epidemic and the subsequent lockdown, i.e., the period since the start of epidemic up to May 7, 2020, to which all interventions outlined in Table 1 applied. As the starting date we arbitrarily used February 1, 2020, in order to cover the real and unknown beginning of the epidemic (Methods; the first three cases apparently became infected and showed symptoms before March 1 and some undetected asymptomatic cases were most likely imported too). The mechanistic character of our model allows for an exact implementation of specific dates of initiation of various interventions and for setting up factors that reflect impacts of these interventions on epidemiological as well as behavioral parameters. We base these factors on data from sociological surveys quantifying behavior before and during the lockdown (Methods, Fig. 2).

As the first step in applying the ABC technique we simulated 200,000 realizations of epidemic dynamics with parameter values generated randomly from prior distributions based on public health data in Czechia and on literature-based data (Fig. 3A; Materials and Methods). The 0.05% of model realizations (100) that best matched the observed age-specific cumulative numbers of confirmed cases were then selected using a rejection-sampling algorithm (Fig. 3C–E). The parameter sets corresponding to the selected realizations form an outer estimate of a posterior distribution of parameter values, the distribution that possesses residual uncertainty in the parameters after the model is confronted with the observed data (Toni et al. 2009; Beaumont 2010; Csilléry et al. 2010) (Fig. 3A). One may think of the selected parameter sets as representing different worlds, the observations of which are compatible with the actual observations, and study the epidemic of Covid-19 in any of these worlds or in all of them simultaneously (Methods).

Fig. 3

Approximate Bayesian computation (ABC). A Prior parameter distributions based on available information are first built for every estimated parameter (blue). Parameter sets randomly generated from these priors are used to run the model; the sets for which the model outcome provides good fit to the observations are selected and used to built a posterior distribution over the estimated parameters (red). C–F Temporal dynamics of the age-specific cumulative numbers of confirmed cases for the baseline scenario, covering the initial phase of epidemic and the following lockdown, for C children (age cohort 0–19 years), D adults (age cohort 20–64 years), and E elderly (age cohort 65+ years), with their sum plotted as panel (B). F Probability density functions for the numbers of confirmed cases by May 7, 2020. In panels (B–E), dark solid curves represent real observations, light (orange in B) curves are results of model simulations run for 100 selected posterior parameter sets, and dark dashed curves are means of those 100 simulations. In panel (F), the dark vertical lines represent the observed numbers of confirmed individuals by May 7, 2020, while the light areas are age-specific density plots for the simulated numbers of confirmed individuals by May 7, 2020, for the 100 selected posterior parameter sets (Color figure online)

3.2 Impact of Timing and Potential Alternative Interventions

Several scenarios alternative to the actually implemented lockdown and motivated by ideas appearing among the general public were considered, run for the selected parameter sets, and compared with the baseline scenario.

The first two scenarios address one of the most important issues in containing an incipient epidemic: the timing of interventions. We find that establishing every single intervention 4 days earlier (scenario R1) or later (scenario R2) than actually happened produces significant differences in the numbers of confirmed cases by May 7, 2020 (Fig. 4). All else being equal, any delay of 4 days in implementing the lockdown thus doubled the number of confirmed cases by May 7, 2020. When deciding on interventions and implementing them time is undoubtedly of essence.

Fig. 4

Lockdown scenario shifted in time. Temporal dynamics of the cumulative numbers of confirmed cases for the lockdown scenario shifted 4 days backward (A and B) or 4 days forward (C and D). Blue solid lines in the left panels represent real observations of the total cumulative numbers, orange lines are the results of model simulations run for the 100 selected posterior parameter sets, and blue dashed curves are means of those 100 simulations. In the right panel, probability density functions for the age-specific numbers of confirmed cases by May 7, 2020, are provided; the dark vertical lines represent the observed numbers of confirmed cases by May 7, 2020, while the light areas are density plots for the simulated numbers of confirmed cases by May 7, 2020, for the 100 selected posterior parameter sets (Color figure online)

Our next two scenarios are motivated by the initial reluctance to order the compulsory wearing of face masks in many countries, and by the existence of groups that deny their utility even today. However, personal protection is not only about face masks. It also comprises increased hygiene and behavioral changes towards protecting oneself. Setting the amount of interpersonal contacts to the sociologically observed levels corresponding to the baseline scenario while neglecting any sort of personal protection ($r_{\mathrm{P}} = 0$), produced many times more confirmed cases by May 7, 2020 (scenario R3, Fig. 5A, B). On the contrary, sticking to the level of personal protection observed in the population yet not limiting interpersonal contacts in any way worsened the epidemic approximately by a factor of three by May 7, 2020 (scenario R4, Fig. 5C, D). The modeled epidemic thus appears less sensitive to changes in contact structure than to changes in the factors reducing infection transmission due to personal protective measures.

Fig. 5

Limiting contacts versus enhancing protection. Dynamics of the cumulative numbers of confirmed cases for A and B the scenario R3 with no use of personal protective measures and C and D the scenario R4 without anyhow limiting contact structure. Legend as in Fig. 4 (Color figure online)

A strategy of sheltering the elderly while letting the remaining population mix relatively freely remains a popular mantra among opinion makers dissenting to population-wide restrictions. Unfortunately, such a strategy does not appear to work (Fig. 6). The explanation for this is, perhaps surprisingly, quite straightforward. The elderly typically represent a relatively small but still sizeable fraction of the population (about 20% of people in Czechia are of age 65 years or older). It is thus impossible to isolate them completely and the extraordinarily high viral load that would occur in the freely mixing younger population would eventually bring the infection even into the protected group. In either of the two scenarios modeling this situation (R5 and R6) the number of affected seniors is (much) higher by May 7, 2020, compared to baseline but even more alarming is the convex rather than concave shape of the curve indicating a continuing rapid increase after that date (Fig. 6).

Fig. 6

Sheltering the elderly does not work. A, C, E Contacts of seniors with the rest of population are reduced by 90%, while contacts within the rest of population are reduced by 10%. B, D, F Contacts of seniors with the rest of population are reduced by 90%, while contacts within the rest of population are reduced by 40%. In both scenarios, personal protection in seniors is set to 90%, while in children and adults only to 50%. Legend as in Fig. 3 (Color figure online)

Finally, an ongoing debate in many countries is whether to leave schools open (or rather whether to open them, under some restrictions or rotation schemes). In a similar vein one may ask whether home office or limiting contacts at work is effective. In our model leaving schools open means on average 4.8 more child-child contacts per day, while no home office keeps on average 5.3 adult-adult work contacts per day (Prem et al. 2017). Leaving all other restrictions as in the baseline scenario, the last two scenarios (R7 and R8) demonstrate the importance of these restrictions. With open schools (Fig. 7A–C) or no home office (Fig. 7B–D) the total number of confirmed cases by May 7, 2020, increases compared to baseline. Since these interventions predominantly affect different parts of the population (children vs. adults), we observe differential effects on the respective age cohorts (Fig. 7B, D). Yet the more interesting observation here is the uncertainty in the qualitative course of infection not observed in any of the previously considered scenarios: whereas some trajectories level off as for the baseline scenario, suggesting a comparatively small effect of the respective interventions, others continue to rise, even though rather linearly than exponentially (Fig. 7A, C).

Fig. 7

Letting schools or work fully open. Dynamics of the cumulative numbers of confirmed cases for A and B the scenario R7 with schools left open and C and D the scenario R8 with not enforcing any home office. Legend as in Fig. 4 (Color figure online)

A way of assessing the effect of (not) implementing an intervention in our multi-world framework is to calculate the impacts of an intervention relative to the baseline scenario when both simulations come from the same parametric world. This way of plotting the results showcases our findings that every 4 days of delaying the interventions result in a doubling of the number of confirmed cases by the end of the lockdown period, that personal protective measures are more effective than just reducing the number of contacts, that sheltering the elderly is not as effective as it may seem, and that leaving schools open or not adopting home office has ambiguous effects (Fig. 8). It also shows that despite different parameter sets that nonetheless produce results matching actual observations, the effects of interventions are comparable across these parametric worlds, providing a robustness check for our model-based assessment of their efficacy.

Fig. 8

Comparing effects of analyzed scenarios. Box plots showing the effects of adopted scenarios in absolute terms (A, C) and relative to the baseline general lockdown scenario (B, D) for the total numbers of confirmed cases expected on May 7, 2020. Codes for different scenarios are explained in the text and Table 5

4 Discussion

Non-pharmaceutical interventions always underlie the first wave of defense against an incipient epidemic. The art of responding to the epidemic is not only in setting up the right restrictions and in the right time but also in ordering them optimally and implementing them as effectively as possible. This all sounds intuitive, and it really is. However, governments still vary in their approaches.

The best scenarios from a theoretical point of view, such as timely and adequate testing and tracing and full tracing coverage, are unfortunately seldom achieved in practice. Mathematical models demonstrate their strength in situations when the optimal strategy is unattainable, showing, for example, what we can expect and how to respond to testing delays of 1 or 2 days or when only 70–80% tracing coverage is doable (Kretzschmar et al. 2020). The major advantage of modeling epidemics is thus in providing quantitative assessments of what kind and/or intensity of interventions is likely (not) enough to prevent serious epidemic impacts, or, more loosely, which strategies of containing epidemic would (not) work.

Two types of models have mostly been used to provide insights into effects of non-pharmaceutical interventions in controlling Covid-19 epidemics. One type, exemplified by Domenico et al. (2020), Ferguson et al. (2020), Kissler et al. (2020), Bubar et al. (2021) and Chang et al. (2021), is a mechanistic prospective model that aims at predicting a future course of epidemic given a set of alternative scenarios. Such models may be calibrated on real observations in a specific country, and typically provide a relative ranking of the considered strategies. Another type of models, exemplified by Brauner et al. (2021), Dehning et al. (2020), Flaxman et al. (2020), Haug et al. (2020) and Liu et al. (2021), are statistical in nature and attempt to reveal effects of interventions by retrospectively analyzing observed time series data on Covid-19 epidemics in many countries at once. Such models typically do not have a means of switching the various interventions on and off, thus losing the ability to explore the impacts of alternative (i.e., not implemented) interaction scenarios. Our model in a sense bridges these two types, providing a mechanistic description of the epidemic with the ability to switch interventions on and off at specific time instants and to modify their intensity while at the same time being calibrated on a robust set of observed data (Davies et al. 2020; Pei et al. 2020).

We used a mathematical model to analyze the first spring 2020 wave of Covid-19 in Czechia. The model allowed us to distinguish between the effects of personal protective measures and interventions reducing social contacts and to set exact dates of implementing such regulations. It structures the population by age (three cohorts: 0–19 years [children], 20–64 years [adults] and 65+ years [seniors]) and type of contacts (four types: at home, at school, at work and other contacts within community). The model was calibrated on age-specific cumulative numbers of confirmed cases, requiring an observation layer that modeled delays from symptoms onset to testing to reporting and a proportion of the population not willing to test themselves despite being symptomatic. We found that (1) delaying the lockdown by 4 days led to twice as many confirmed cases by the end of the lockdown period, (2) personal protective measures such as wearing face masks and using hand sanitizers were more effective than only reducing social contacts, (3) sheltering only the elderly while letting the remaining population live relatively normally was not a viable strategy, and (4) leaving schools open and not adopting home office appeared to be risky measures.

A key and unique feature of our model is its use of sociological data quantifying the degree of compliance with the interventions, that is, the degree of interpersonal contact limitation in various environments as well as the degree of compliance with using personal protective measures. Therefore, we need not speculate on the degree of compliance in the society but rather work with the actual level of the reduction in social contacts relative to the preceding non-epidemic state. These proportional reductions in contacts enter our model as fundamental system characteristics and serve as benchmarks, towards which we test alternative intervention strategies.

Another unique feature of our modeling approach is how we deal with uncertainty. Although our model is deterministic we calibrate it via a stochastic technique, known as Approximate Bayesian Computation (ABC) (Toni et al. 2009; Beaumont 2010; Csilléry et al. 2010). In this way we generate a number of posterior parameter sets (or different worlds) that all produce outcomes matching the actual observations of the age-specific numbers of confirmed cases, and examine how the epidemic behaves in those worlds when the interventions change. We emphasize that even though our model is calibrated for the epidemic in the Czech Republic it can straightforwardly be recast for any other country or region if relevant data on behavioral and population characteristics are available.

Despite the obvious fact that any delay in implementing interventions aimed at slowing down the epidemic has a negative impact, many politicians and public health officers tend to downplay its importance. Here we quantified this effect for the lockdown applied during the first spring 2020 wave of Covid-19 in Czechia, showing that a delay of 4 days doubles the number of confirmed cases by the end of lockdown period. All else being equal, this also means doubling the number of deaths during that period, as the change affects all age cohorts equally. Clearly, this exact quantitative relationship cannot be expected to hold for any epidemic wave and in any country but it shows that the effect is far from negligible and has to be taken seriously. In light of this result it is quite concerning to see the specifics of lockdown implementation in Czechia: it has been sadly common that interventions announced as a response to a worsening epidemic situation were actually implemented with a delay of 1 or 2 weeks. A similar exploration of Covid-19 epidemic in the United States between March 15 and May 3, 2020, suggested that implementing their interventions 1 or 2 weeks earlier would have reduced the numbers of confirmed cases and deaths by 50% and 90%, respectively (Pei et al. 2020). This further underscores the need to act as quickly as possible.

Any epidemic can only be brought under control by interrupting the chain of transmissions. This can be done by limiting social contacts and/or by adopting personal protective measures. Since locking ourselves at home or wearing full protective suits for 2 or 3 weeks are obviously not viable options, it is important to examine the relative effect of limiting contacts versus enhancing protection. In this light it is surprising how long people (and even some governments) resisted using face masks in many parts of the world, with Sweden recommending wearing face masks only in December 2020 (BBC 2020). On the other hand, countries such as Italy or Austria insisted on wearing face masks also during the relatively “peaceful” summer months. A positive impact of wearing face masks on limiting spread of Covid-19 has been demonstrated by a number of modeling studies (Eikenberry et al. 2020; Stutt et al. 2020). We add to this knowledge through finding that adhering to personal protective measures is likely more effective than just reducing social contacts. It is then not surprising that our simulations also suggested that switching the dates of implementing wearing face masks and closing schools in Table 1 resulted in only half the amount of confirmed cases by May 7, 2020 (results not shown). For any respiratory epidemic to come in the future using face masks and other personal protective measure should therefore be among the first interventions implemented.

The idea that sheltering of elderly, or more broadly any high-risk portion of the population, while letting the rest live more or less normally, is sufficient to prevent deaths still has its followers. We show here that it is deeply flawed. First, in Czechia and many other European countries the age cohort of age 65 years or more forms at about 20% of the population, of which only a small part lives in senior houses while many others live in multi-generational families, so it is difficult to shelter them successfully. Second, leaving the rest of the population behaviorally unrestrained would soon result in a high virus prevalence in this group. The virus would then sooner or later percolate into the elderly group despite its relatively high protection. Our simulations show that even a 40% reduction in social contacts among the rest of the population would not help in preventing the infection of the elderly. Moreover, instead of the decelerating by May 7, 2020, the scenario of sheltering the elderly exhibits a continuing exponential growth. This scenario is thus not a viable option, even more so if we take into consideration the potential impact of such isolation on the mental state of the elderly: Krendl and Perry (2021) found that older individuals felt higher depression and greater loneliness during pandemic.

Finally, an issue commonly discussed in both country-specific studies (Davies et al. 2020) and studies spanning more countries (Brauner et al. 2021; Flaxman et al. 2020), is whether to leave schools open or to close them. The closing of schools was among the first interventions in many countries (Brauner et al. 2021; Flaxman et al. 2020), presumably because of existing pandemic plans tailored to influenza. We considered scenarios in which schools remained open or alternatively home office was absent, with all other interventions kept as during the lockdown, with ambiguous results. Whereas in many parametric worlds an effect of leaving schools open or absenting home office appears negligible, in many remaining worlds an increase in the number of confirmed cases becomes linear rather than decelerating by May 7, 2020. This is worrying since a significant increase in the prevalence among younger age cohorts sooner or later manifests itself in the older ones, inevitably resulting in more deaths. Other studies are likewise ambiguous on this issue, and can be generally divided into those that rank closing schools among the most important interventions (Brauner et al. 2021, and references therein) and those that claim a relatively small effect of closing schools (Davies et al. 2020; Flaxman et al. 2020; Iwata et al. 2020; Viner et al. 2020). Discordant results of these studies may be attributed to using different type of models, which calls for caution when reading and interpreting such studies, including our one. For example, both Brauner et al. (2021) and Flaxman et al. (2020) admit that some of their conclusions might have been affected by the fact that the analyzed interventions were applied in a specific order or close to each other.

There are likely spin-off effects of various measures, such as when schools are closed kids spend more time with grandparents. After closing schools in Czechia, many students enjoyed meeting themselves in shopping malls. The answer to this was quick: a few days later wifi had to be shut down in shopping malls and all resting places and food-courts were closed. Effects of closing schools or more frequent home office on home contacts are more questionable, since contacts in families are commonly long enough and close enough for transmission also under common conditions. In fact, in the Czech Republic, the lockdown also banned neighbor contacts or contacts within wider families (family visits). In summary, it is difficult to assume real extent of spin-off effects.

Although the majority of compartmental epidemic models, including those used to describe Covid-19, assume linear rates of leaving specific model classes (thus corresponding to exponentially distributed residence times), others are formulated as delay differential equations and consider the residence times as fixed (Smith 2010). For example, acknowledging that an epidemic cannot persist if most infected hosts are isolated in a sufficiently short time, Ruschel et al. (2019) and Young et al. (2019) developed an SIQ model with an isolated class and derived for it a critical proportion of infected individuals and a critical time within which an isolation has to be made in order to stop an epidemic. Moreover, whereas our model is discrete in time (hence involves a time step) and uses linear rates of leaving specific model classes, the model of Young et al. (2019) is continuous in time and combined linear (recovery and cessation of immunity) and fixed (latent and isolation periods) delays. Disregarding that in the latter model recovered individuals will become susceptible again, it would be worth comparing various approaches to modeling time delays using a common data set on a real epidemic. Reality is arguably somewhere in between, and no golden rule of modeling time delays appears to exist.

Here, we use a deterministic epidemic model, complemented by the stochastic Approximate Bayesian Computation technique that allows one to obtain several alternative parametric sets with all of which the model reasonably fits the observed epidemic course. The use of these parametric sets for running the model under-examined scenarios thus introduces a kind of variability imposed by parametric uncertainty. Alternative approaches considering stochastic models with perturbed parameter values and/or providing individual parameter estimates are frequently used (Kucharski et al. 2020; Rozhnova et al. 2021). Again, a methodological study quantifying contribution of each of these (and also other) sources of uncertainty to the overall model sensitivity in the context of a real epidemic is naturally warranted.

Lockdown is an expensive measure and should generally be applied only when milder interventions fail. With a limited knowledge of Covid-19 and the resulting fear during the first spring 2020 wave, many countries adopted a form of lockdown, allowing us to take it as a baseline scenario with which to compare alternative strategies and thus help assess interventions for the other epidemic waves. We used a well-defined model that allowed us to play with the adopted interventions and at the same time calibrate the model on the lockdown situation. Our results suggest that a combination of a timely application of personal protective measures, a limitation of contacts and an effective testing (and contact tracing) constitute a sound strategy. This may sound trivial but mathematical models can tell us how much is not enough, what does not work and what appears to be a necessary minimum to keep the society alive and to protect its health.