On the Parametrization of Epidemiologic Models—Lessons from Modelling COVID-19 Epidemic

Numerous prediction models of SARS-CoV-2 pandemic were proposed in the past. Unknown parameters of these models are often estimated based on observational data. However, lag in case-reporting, changing testing policy or incompleteness of data lead to biased estimates. Moreover, parametrization is time-dependent due to changing age-structures, emerging virus variants, non-pharmaceutical interventions, and vaccination programs. To cover these aspects, we propose a principled approach to parametrize a SIR-type epidemiologic model by embedding it as a hidden layer into an input-output non-linear dynamical system (IO-NLDS). Observable data are coupled to hidden states of the model by appropriate data models considering possible biases of the data. This includes data issues such as known delays or biases in reporting. We estimate model parameters including their time-dependence by a Bayesian knowledge synthesis process considering parameter ranges derived from external studies as prior information. We applied this approach on a specific SIR-type model and data of Germany and Saxony demonstrating good prediction performances. Our approach can estimate and compare the relative effectiveness of non-pharmaceutical interventions and provide scenarios of the future course of the epidemic under specified conditions. It can be translated to other data sets, i.e., other countries and other SIR-type models.


Introduction
Predicting the spread of an infectious disease is a pressing need as demonstrated for the present SARS-CoV-2 pandemic. Due to the worldwide high disease burden, a plethora of mathematical epidemiologic models was proposed. This includes auto-regressive time series methods, Bayesian techniques, and application of deep learning methods, but also mechanistic models and hybrid models combining some of these approaches (see [1] for a review). Among mechanistic models, based either on agents [2][3][4] or on compartments [5], the most commonly proposed and published model type is the classical SIR (S = susceptible, I = infected, R = recovered) type compartment model, which was presented with different modifications often considering further aspects and details such as disease states, age structure, contact patterns, and intervention effects. Major aims of these models are to predict (1) the dynamics of infected subjects, (2) requirements of medical resources during the course of the epidemic, or (3) the effectiveness of non-pharmaceutical intervention programs (NPI) [6,7]. Examples for models addressing these aims are a SECIR model (E = exposed, C = cases) proposed by Khailaie et al. [8], a SEIR type model proposed by Barbarossa et al. [9], models from The Robert-Koch institute [10], and from Dehning et al. [11].
A good prediction performance does not only depend on the precise structure of the model but on its parametrization. This, however, is a non-trivial and often underestimated task due to the following issues applicable to other infectious diseases as well: Viruses 2022, 14, 1468 2 of 37 key epidemiologic parameters are often unknown or only known within ranges. Therefore, parametrization based on observational data is a common approach. However, reported official data bases are heterogeneous and often biased due to (1) lag in reporting of cases/events [12], (2) changing testing policy either due to limited testing capacities, which might depend on the pandemic situation itself or by changing risk profiles of people to be tested (e.g., defined risk groups, dependence on symptoms, degree of prophylactic testing) [13], and (3) incompleteness of data [14]. Moreover, parametrization depends on further epidemiologic issues to be considered, comprising (1) changing age-structure of the infected population with impact on symptomatology, hospital or intensive care requirements and mortality, (2) spatial heterogeneity of the spread of the disease driven by local conditions and outbreaks, (3) new pathogen variants becoming prevalent, (4) nonpharmaceutical interventions continuously updated in response to the pandemic situation, and finally, (5) the progress of vaccination programs and its effectiveness.
Due to these interacting complexities, it is close to impossible to construct a fully mechanistic model covering all these aspects in parallel. Therefore, we here propose a framework of epidemiologic model parametrization, which accounts for these issues in a more phenomenological, data-driven way applicable even for limited or biased data resources.
In detail, we here propose to integrate a mechanistic epidemiologic model as a hidden layer into an input-output non-linear dynamical system (IO-NLDS), i.e., the true epidemiologic dynamics cannot be directly observed. This allows distinguishing between features explicitly modelled (in our case different virus variants, vaccination) and changing factors of the system which are difficult to model mechanistically (in our case changes of contacts, e.g., due to non-pharmaceutical interventions or changing contact behavior, changing age-structure of the infected population and changes in testing policy, in the following abbreviated as NPI/behavior). These factors are imposed as external inputs of the system. We then estimate parameters by a knowledge synthesis process considering prior information of parameter ranges derived from different external studies and other available data resources such as public data. We are thus going beyond previous modeling approaches that only used point estimates for parameters [15,16]. Specifically, we use Bayesian inference for the parameter estimation, which could also be time-dependent. We analyze the structure of available public data in detail and translate it to model outputs linked by an appropriate data model to the hidden states of the IO-NLDS, i.e., the epidemiologic model. We demonstrate this approach on an example epidemiologic model of SECIR type for SARS-CoV-2 and data of Germany and Saxony, but our method can be translated to other countries, other models or even other infectious disease contexts.

General Approach
We consider input-output non-linear dynamical systems (IO-NLDS) originally proposed as time-discrete alternatives to physiological pharmacokinetic and -dynamics differential equations models [17]. This class of models couples a set of input parameters such as external influences and factors with a set of output parameters, i.e., observations by a hidden model structure to be learned (named core model in the following). This coupling is not necessarily fully deterministic, i.e., data are not required to represent directly state parameters of the model. This represents a major feature of our approach in order to account for different types of biases in available observational time series data.
We here demonstrate our approach by using an epidemiological model as core of the IO-NLDS. Non-pharmaceutical interventions, changes of testing policy, age distributions and severity of disease courses were phenomenologically modelled by external control parameters imposed on the epidemiologic model via the input layer of the IO-NLDS. Random influxes of infected subjects e.g., by travelling activities or outbreaks are also considered by this approach. Number of reported infections, intensive care (IC) cases, and deaths are considered as output parameters not directly representing the hidden states of Viruses 2022, 14,1468 3 of 37 the model due to several data issues including reporting delays. The model is then fitted to data by a full information approach, i.e., all data points were evaluated by a suitable likelihood function.
The single steps of this process are explained in detail below.
Assumptions of the core model We adapted a standard SECIR model (SECIR = susceptible, exposed, cases, infectious, recovered) for pandemic spread. We introduced an asymptomatic compartment in order to account for infected patients, which do not have symptoms, a common condition of SARS-CoV-2 infection. A compartment of patients requiring intensive care (IC) was added to model respective requirements and we distinguished between deceased and recovered patients.
We subdivided most of the compartments into three sub-compartments with first order transitions to model time delays. Transition rates between sub-compartments are the same for each respective compartments for the sake of parsimony. This approach is extensively used in pharmacological models [10]; it is equivalent to a Gamma-distributed transit time [11]. To allow for two concurrent virus variants with differing properties, compartments of asymptomatic and symptomatic infected subjects are duplicated. This allows us, for example, to simulate the take-over of the more infectious B.1.1.7 (Alpha), and later, B1.617.2 (Delta) variant observed e.g., in all European countries [12].
The general scheme of the IO-NLDS system is shown in Figure 1. We make the following assumptions: 1.
The input layer consists of external modifiers influencing (1) reporting policy (e.g., changing testing policy), (2) rates of infections (affected by non-pharmaceutical interventions, age structure, influx of cases), and (3) risks of severe disease conditions such as IC requirements and deaths, also depending on the changing age structure of infected subjects.

2.
The output layer of observable data is linked to the hidden layers of the core model by specific data models (see later).

3.
Susceptible, non-infected people (Sc): We assume that 100% of the population is susceptible to infection at the beginning of the epidemic. 4.
The latent state E comprises infected but non-infectious people. 5.
The asymptomatic infected state IA has three sub-compartments (I_(A,1), I_(A,2) and I_(A,3)). From I_(A,1), transitions to the symptomatic state or the second asymptomatic state are possible. From I_(A,2), only transitions to I_(A,3) and then to the recovered state R are assumed. 6.
Both cases I_A and I_S contribute to new infections but with different rates to account for differences in infectivity and quarantine probabilities. 8.
The compartment C represents critical disease states requiring intensive care. We assume that these patients are not infectious due to isolation. Again, the compartment is divided into three sub-compartments, C_1, C_2, and C_3. In C_1, a patient can either die or transit to C_2, C_3, and finally, R. 9.
Patients on the recovered stage R are assumed to be immune against re-infections. 10. We duplicate the compartments E, I_(A,1), . . . , I_(A,3), I_(S,1), . . . , I_(S,3) to account for two concurrent virus variants. We assume different infectivities for the two variants. All other parameters are assumed equal. No co-infections are assumed. The input layer consists of external modifiers including parameter changes due to changes in testing policy, non-pharamceutical interventions, and age-structures. The output layer is derived from respective hidden layers via stochastic relationships (see later). The output layer is compared with real-world data. The superscript Mu denotes new virus variants.
These assumptions are translated into a difference equation system (see Appendix A). Model compartments and their properties are explained in Table 1. Symptomatic infected state 1. Can either become critical, i.e., transits to C 1 with probability p crit and rate r 6 or stays sub-critical with probability 1 − p crit and transits to I S,2 with rate r 5 All model parameters of the model are described in Table 2. Complete dynamics of the epidemic in Germany is shown in Figures 2 and 3.

Input Layer
The input layer represents external factors acting at the SECIR model, effectively changing its parameters [42]. We define step functions b 1 and b 2 as time-dependent input parameters modifying the rate of infections caused by asymptomatic, respectively symptomatic subjects. To identify dates of change, we used a data-driven approach based on a Bayesian Information Criterion informed by changes in non-pharmaceutical interventions for Germany based on Government decisions, changing testing policies as well as events with impact on epidemiological dynamics such as holidays or sudden outbreaks. Details can be found in Appendix B.
We also accounted for changes in the probabilities of critical courses and mortality, which can be explained by changes in testing policies covering asymptomatic cases to a different extent (for example symptomatic testing only vs. introduction of screening tests, e.g., rapid antigen tests), respectively shifts in the age-distribution of patients or changes in patient care efficacy (new pharmaceutical treatment, overstretched medical resources). Again, this is implemented by step functions p crit , respectively p death . Number of steps are determined on the basis of a Bayesian Information Criterion. Details can be found in the Appendix B as well as in Table A5 from Appendix I and Table A8 from Appendix J. The parameter P S,M represents the percentage of reported infected symptomatic subjects in relation to all symptomatic subjects. This value is assumed to be constant (50%) in the present version of the model. We describe the parameters defining the input layer in Table 3. Table 3. Parameters used to define the input layer. These parameters were used to empirically model changing NPIs or changing contact behavior, changes in testing policies and changing age-structures during the course of the epidemic. Respective input functions constitute the input layer of our IO-NLDS model.  However, also for the considered time intervals several sources of bias need to be considered. We handled these issues as explained in the following:

Parameter
Infected cases: We first smoothed reported numbers of infections with a sliding window of seven days centered on the time point of interest to control for strong weekly periodicity. We assume that these numbers correspond to a certain percentage P S,M of symptomatic patients. This is justified by the fact that the majority of reported infected cases develop symptoms (about 85% according to the RKI [43]), but there is also a large amount of asymptomatic cases (approximately 55-85% of infections [37][38][39]. In the present model, we assume P S,M as constant. The exact equation linking states of the SECIR model with the measured numbers of infected subjects I S,M can be found in Appendices B and C Equation (A7). We further accounted for delays in the reporting of case numbers by introducing a log-normally distributed delay time as explained in Appendix C.
Critical cases: The number of critical COVID-19 cases (DIVI reported ICU) is available since end of March 2020 [44]. We assumed that these data are complete since 16 April 2020 when reporting became mandatory by law in Germany. Earlier data were up-scaled from the number of reporting hospitals to the number of ICU-beds of all hospitals according to the reported ICU capacity available for 2018. We coupled the sum of the critical subcompartments C i (i = 1,2,3) to these numbers directly.
Deaths: Deaths are reported by the RKI but daily reports do not reflect true death dates, which needs to be accounted for. Available daily death data of the RKI are retrospectively updated with a delay between true death date and reported date (death report delay-DRD). We assume that the DRD is normally distributed with an average of 7.14 days and a standard deviation of 4 days as reported by Delgado et al. [45]. Details can be found in Appendix C.
Occurrence of B.1.1.7 variant: In January 2021, the variant B.1.1.7 became endemic in Germany and quickly replaced all other variants. Onset of this development was modeled by an instantaneous influx of 5% of newly infected subjects into the E Mu compartment on 26 January 2021 estimated from published data [46].

Parametrization
We carefully searched the literature to establish ranges for our model parameters. These ranges are used as prior constraints during parametrization of our model ( Table 2). Justification of prior values is provided in Appendix H. Parameters are then derived by fitting the predictions of the model to reported data of infected subjects, ICU occupation, and deaths using the link functions of model and data explained in the previous section. This is achieved via likelihood optimization. Likelihood is constructed using the same principles as reported [47]. In short, the likelihood consists of three major parts, namely the likelihood of deviations from prior values, the likelihood of residual deviations from the data, and a penalty term to ensure that model parameters are within the prescribed ranges, as explained in Appendix D. We follow a full-information approach intended to use all data collected during the epidemic as explained in Appendices E-G. As a result, our model fits well to the complete dynamics of the epidemic in Germany in the above mentioned time period (Figures 2 and 3).
To ensure identifiability of parameters, we checked a number of parsimony assumptions. For example, we assumed that the dynamical infection intensities of asymptomatic (b 1 ·r 1 (t)) and symptomatic subjects (b 2 ·r 2 (t)) are proportional with factor rb 1,2 . We also determined Bayesian Information criteria (BIC) for different partitioning numbers of the external jump functions (N crit and N death ) to keep these as small as possible. Details can be found in Appendix B.
Likelihood optimization is achieved using a stochastic approximation of an estimationmaximization algorithm (SAEM) [48]. The algorithm is based on a stochastic integration of marginal probabilities without using likelihood approximations such as linearization or quadrature approximation or sigma-point filtering [17]. Confidence intervals of model predictions are derived by Markov-Chain-Monte-Carlo simulations, i.e., alternative parameter settings were sampled from the parameter space around the optimal solution (Appendices B, F and G). We use these parameter sets to simulate alternative epidemic dynamics. This resulted in a distribution of model predictions from which empirical confidence intervals are derived.

Implementation
The model and respective parameter estimations are implemented in the statistical software package R from which external publicly available functions are called. The model's equation solver is implemented as C++ routine and called from R code using the Rcpp package. The code and data for simulation of the output layers using the reported parameter settings will be made available via our Leipzig Health Atlas: (https://www. health-atlas.de/models/40, accessed on 26 June 2022) and GitHub (https://github.com/ GenStatLeipzig/LeipzigIMISE-SECIR, accessed on 26 June 2022) [49].

Explanation of Epidemiologic Dynamics
We used the full data set to explain the course of infections, ICU occupations, and deaths between 4 March 2020 and 29 March 2021 in Germany. A total of three parameters were assumed time-dependently, namely Infectivity b 1 and the probability of a critical disease course (p crit ) and death (p death ). We identified nine fixed and 19 empirically identified time points of NPI/behavioral changes (Table A1 from Appendix I). Regarding p crit and p death , we identified 18 respectively 19 time steps (See Tables A2 and A3 from Appendix I and Table A7 from Appendix J). Throughout the epidemic, we observed a good agreement of our model and incident ( Figure 2) and cumulative data ( Figure 3). Corresponding residual errors are provided for all observables (Table A4 from Appendix I). As shown in Table A2 from Appendix I, we estimated 14 static and three dynamically changing parameters using 1170 data points (390 daily measurements of registered cases, registered deaths and ICU occupancy) for Saxony as well as for Germany. Comparison of IO-NLDS model (magenta curve) and data (thin grey curves = raw data, solid black curve = data averaged by sliding window) is provided in the upper column. A good agreement is observed (shaded area = prediction uncertainty, vertical lines = changes in NPI/contact behavior). The middle row represents the corresponding input layer, i.e., the estimated time course of the time-dependent input parameters, namely infectivity and probabilities of critical disease course and death. Time steps correspond to the lines of changing NPI/contact behavior as displayed in the upper row. In the lower row, we present percentages of B.1.1.7 among infected subjects (first figure), subjects older than 80 years among infected corresponding to high death tolls (second), and subjects in the age categories 35-59, respectively 60-79 among critical cases (last figure of last row).

Parameter Estimates and Identifiability
Parameter estimates of the SECIR model are presented in Table 2, while those required to define the input layer are presented in Table 3 and Tables A1 and A3 from Appendix I. For those parameters for which we used prior information for fitting purposes, we compared the respective expected posteriors with their best priors (see Figure 4). Statistics are provided in Table A5 from Appendix I. No significant deviations between expected values of posteriors and priors were detected. All relative errors of parameters of the SECIR model are smaller than 10% demonstrating excellent identifiability of all epidemiologic parameters. As expected, identifiability of the external control functions is reduced. Largest standard errors of steps are in the order of 70% still demonstrating reasonable identifiability (Table A3 from Appendix I).   Table A8 from Appendix J.

Plausibilization of Estimated Step Functions of Infectivity
We estimated the infectivity as an empirical step function through the course of the epidemic. This step function should also roughly reflect NPI effectivity. We therefore compared our infectivity step function with the Governmental stringency index of NPI as estimated on the basis of Hale et al. [50]. Results are displayed at Figure 5 and revealed a reasonable agreement.

Model Predictions
We regularly used our model to make predictions regarding the future course of the epidemic. Predictions were specifically made for the Free State of Saxony, a federal state of Germany and were published at the Leipzig Health Atlas [49]. We here present comparisons of our predictions with the actual course of the epidemic for two scenarios to demonstrate utility of our approach. Parameter values for Saxony were obtained in the same way as for Germany restricting available data of infected subjects, ICU cases, and deaths to this state. Estimated parameter values are presented in Tables A6-A8 from Appendix J. While Saxony was almost spared from the first wave of SARS-CoV-2 in Germany, the second wave hit the country particularly hard resulting in the highest relative death toll of all German states (1 out of 400 inhabitants of Saxony died from COVID-19 during the second and the immediately following B.1.1.7-driven third wave). The second wave was on its peak in the middle of December 2020. A hard lock-down was initiated at this time including closure of schools, prohibition of all team-based leisure activities, and night-time curfew. We were asked by the government to estimate the length of lock-down required to break the second wave. Stringency of lock-down was comparable to the first wave. Thus, we simulated four scenarios: an optimistic assumption of a lock-down efficacy of 60% reduction in infectivity, a more realistic scenario with 40% reduction, a pessimistic assumption of only 20% lock-down efficacy, and finally, 0% reduction (no lock-down) as control scenario. Results are shown in Figure 6 and revealed a good agreement of our prediction with the actual course for the 40% scenario considered likely.
At the beginning of February 2021, the second wave was broken in Saxony and first relaxations of NPIs were conducted. At this time, the more virulent B.1.1.7 variant became endemic in Germany. At 14 February, the true percentage of the B.1.1.7 variant was unknown due to lack of sequencing capacities. Moreover, there were uncertainties with respect to the increase in infectivity by the B.1.1.7 variant. We therefore simulated three scenarios (optimistic: 10% initial proportion of B.1.1.7, infectivity increased by factor 1.7, expected: 20% initial proportion, 1.8-times increase in infectivity, pessimistic: 30% initial proportion, 2-times increase in infectivity). Results are shown in Figure 7. The actual course of the epidemic was close to the pessimistic scenario, i.e., the second wave was directly followed by a third wave due to the B.1.1.7 variant. Indeed, later data revealed that the proportion of B.1.1.7 was already close to 30% (pessimistic assumption) at the time the simulation was performed. Moreover, our model correctly predicted the variant replacement by B.1.1.7. Figure 6. Comparison of predicted and observed decline of the second wave in Saxony according to the initiated lock-down. Our model was used to fit the observed data until 21 December 2020 (shown as grey curve (raw data) and black curve (smoothed) of reported test-positives). Estimated step functions b 1 and b 2 describing the infectivity of asymptomatic and symptomatic subjects were reduced by 0% (yellow: no lock-down = control scenario), 20% (green: pessimistic scenario), 40% (blue: realistic scenario), and 60% (magenta: optimistic scenario) to simulate four scenarios of the future course of the epidemic under lock-down conditions. The observed numbers of test-positives after the 21 December 2020 are shown in red (light red = raw data, dark red = smoothed) closely followed the expected scenario of 40% lock-down efficacy. Shaded areas represent 95% prediction intervals. The predictions and parameters were reported in our regular bulletin deposited at Leipzig Health Atlas, ID: 85AH9JMUFM-4. black curve = smoothed). Three scenarios were simulated differing in assumed initial proportion of B.1.1.7 which was not exactly known at this time point (10%, 20%, and 30%, respectively) and in the assumptions regarding increased virulence of B.1.1.7 (parameter mur = 1.7, 1.8 and 2, respectively). Predicted course of subjects infected with the respective variants are shown as shaded areas. The observed total numbers of testpositives (light red = raw data of reported testpositives, dark red = smoothed curve) closely followed the pessimistic scenario 3. Lower row: When comparing the proportion of B.1.1.7 as retrieved from [52] from 18 July 2021, initial proportion of B.1.1.7 was indeed close to that assumed for scenario 3. Blue curves represent 95% confidence intervals of the B.1.1.7 proportion predicted for the different scenarios. All predictions were reported in our bulletin at the 20 February 2021 deposited at our Leipzig Health Atlas [53].

Discussion
In this paper, we propose a method of parametrization of COVID-19 epidemiologic models and applied it to an extended SECIR-type model to explain the course of the epidemic in Germany and one of its federal state, the Free State of Saxony. Moreover, we demonstrated how the model can be used to make relevant predictions, which could be validated on the basis of subsequent observational data.
A key idea of our approach is the embedding of differential equations-based epidemic modelling into an input-output dynamical system (IO-NLDs). This has two major advantages. First, the approach allows combining explicit mechanistic models of epidemic spread and phenomenological considerations of external impacts on model parameters via the input layer. This allows parametrizing models of different complexity. For example, in our model we non-explicitly considered the effect of age structures of the diseased population by time-dependent input parameters such as probabilities of critical disease courses and deaths. This could easily be replaced for example by age-structured models. We believe that such a combined empirical/mechanistic approach is well suitable to address the complexity of COVID-19 epidemic dynamics for which it is impossible to consider all relevant mechanisms explicitly and in parallel.
The second major advantage of our approach is that we assumed a non-direct link between state parameters of the embedded SECIR model and observables. This allows interposing a data model considering known biases of the available data resources. We aimed at identifying relevant bias sources as far as possible and considered them in our proposed data models. However, these data models could be subjected to changes in the future for example if better data of COVID-19-related death will be released. Improved data models could be easily integrated into our framework.
Note that the IO-NLDS implementation translates the embedded differential equations model to a discrete scale (i.e., days in our case), which however appears to be sufficient for describing an epidemic.
We also want to note that the SECIR model used here is by far neither unique nor the most comprehensive one. For example, The Robert-Koch institute developed a model for the purpose of estimating the effect of different vaccination strategies which could easily be included into our SECIR-type models [54]. Although integration of differential equations-based models into our IO-NLDS context is more straightforward, our approach is also applicable to agent-based models. In general, the aspect of parameter estimation of such models is underdeveloped in view of the highly biased data resources used and to our knowledge, no generic concept was proposed so far.
Based on our IO-NLDS formulation and data models, we parametrized our model on the basis of data of infection numbers, critical cases, and deaths available for Germany and Saxony. Here, we chose a full-information approach considering all data in between start of the epidemic 4 March 2020 to 29 March 2021. We also applied a Bayesian learning process by considering other studies to inform model parameter's settings. Thus, we combine mechanistic model assumptions with results from other studies and observational data. This approach is very popular in pharmacology [55] but despite its importance it is yet rarely applied in epidemiology [11]. It resulted in a complex likelihood function, which is optimized on the basis of Markov-Chain Monte Carlo (MCMC) algorithms, as we described in Appendices B and F. If the likelihood has a unique maximum, most of the samples eventually accumulate in its vicinity after a certain number of "burn-in" steps. This allows an effective MCMC search of the best parameter estimates as well as approximations of their standard errors (standard deviations of the sample) and the degree of overfitting. However, if parameters are interdependent, MCMC algorithm samples manifolds of alternative solutions, resulting in large standard errors of the overfitted parameters. We successfully addressed this issue by a modified version of Maire's algorithm [56]. Central to this approach is the idea that the proposal distribution adapts to the target by locally adding a mixture component when the discrepancy between the proposal mixture and the target is deemed to be too large. In other words, this algorithm samples multidimensional parameters sets, approximating it as a mixture of multivariate Gaussian distribution. Such approaches enable adequate sampling of model parameters and detection of overfitting as well as of multiple local maxima of the likelihood. Our results revealed small standard errors indicating lack of overfitting, see Tables A1 and A3 from Appendix I, Tables A6  and A7 from Appendix J. We also applied rigorous information criteria to limit the number of steps of our input functions. As a consequence, it was possible to identify both, the fixed parameters of the SECIR model and the time-variable input functions representing changing NPI/contact behavior and age-structures.
Model parametrization resulted in a good and unbiased fit of data for the period considered for Germany. Fixed parameter values of the SECIR model did not significantly deviated from their prior values if available. It required 18 respectively 19 steps of changes of the probabilities to develop critical stage and to die respectively. A total of 13 intensification and 15 relaxation events were necessary to describe the epidemic dynamics over the time course of observations. Estimated infectivity roughly correlated with the Governmental Stringency Index [51]. We regularly contributed forecasts of our model to the German forecast Hub [57].
We also demonstrated utility of our model by several mid-term simulations of scenarios of epidemic development in Saxony, a federal state of Germany. We could show that predictions of reported infections were in the range of later observations for scenarios considered likely.
As future extensions and improvements of our model, we will consider stochastic effects on a daily scale, for example to model random influxes of cases or to model random extinctions of infection chains. These effects are relevant to be considered in times of low incidence numbers such as those observed in Germany in the summers 2020 and 2021. Our IO-NLDS framework is well suited to implement such extensions [17].
In future versions of our model, we will also include age-structures and implement a vaccination and waning model in analogy to other research groups. In the current version of the model, we assumed a constant proportion of symptomatic patients reported as infected. This does not consider for example changing testing policies (i.e., symptomatic vs. prophylactic testing). We plan to refine our model in this regard in the future. Finally, we will consider the Delta and Omicron variants emerging in 2021 [53] in the next update of our SECIR model.
In summary, the primary focus of the paper is an adequate parametrization of epidemiological models on the basis of complex, possibly biased data, as well as its coupling with structurally unknown dynamical external influences. This approach allows for a clear separation of mechanistic model compartments from random or time-dependent non-mechanistic influences and biases in the data. We believe that this approach is useful not only for the parametrization of the SECIR model presented here but also for other epidemiologic models including other disease contexts and data structures.  Institutional Review Board Statement: Ethical review and approval were waived for this study as only published data from official sources was used (see manuscript for detailed description of data sources).
Informed Consent Statement: Patient consent was waived for this study as only published data from official sources was used (see manuscript for detailed description of data sources).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Equations of the SECIR Model
We here present the equations of the SECIR model serving as hidden layer of our input-output non-linear dynamical system (IO-NLDS). To fit in this context, the ordinary differential equations of the SECIR model are approximated by a difference equation system describing compartment changes at single days, i.e., time-steps ∆t equals one day. Compartments of the model are explained in Table 1 of the main paper. Parameters are explained in Table 2 of the main paper.
The functions X represent decisions regarding the further disease course defined below.
Analogously, the equations for the concurrent variant compartments are as follows: In some of the compartments, there is a decision between recovery or deterioration of the disease course. For example, patients in C 1 can either die or recover with different rates r i and r j . To model this decision, we split the number of patients in the compartment by a respective probability p and determine the decision factor X as follows.
To start the epidemic in Germany, we assumed an entry of infected cases by a linearly decreasing function starting at 4 March 2020 and becoming zero at 10  at 26 February 2021 was fitted to be 5.3% of the corresponding sum of the normal compartments at the same day.

Appendix B. Input Layer
The input layer of our IO-NLDS is designed to describe the effects of non-pharmaceutical interventions (NPI) and other impacts on infectivity such as behavioral changes, changes in age-structure, testing policy, seasonal effects, or larger outbreaks (abbreviated as NPI/ contact behaviour). Since these effects typically affect contact matrices in different ways, we model this phenomenologically by time-dependent reductions or increases of infection rates caused by symptomatic and asymptomatic subjects as explained in Equation (A1).
We make the following assumptions for non-pharmaceutical interventions: 1. We introduce the relative infectivity function b(t), which changes according to NPI/contact behavior modifications. This is modelled by a linear increase (in case of relaxation) or decrease (in case of tightening) within a fixed time Del tr of two days. Otherwise, b(t) is constant. We denote {T tr,s } N tr s=1 as the time points with changes in non-pharmaceutical interventions with N tr the total number of time points with changes. We collected dates of changing non-pharmaceutical intervention measures for Germany based on government decisions, changing testing policies as well as events with impact on epidemiological dynamics such as holidays and sudden outbreaks(such as thin peak of new infections in June affected mostly workers of the meat industry). We also assumed additional time points with changes determined by BIC.
Again, for the sake of parsimony, we assume that the relative infection intensities of asymptomatic (b 1 (t)) and symptomatic subjects (b 2 (t)) are the same, hence, respective proportionality rb 1,2 is constant and estimated during model fitting.

Thus, b(t) is defined as follows:
maceutical interventions with the total number of time points with changes. We collected dates of changing non-pharmaceutical intervention measures for Germany based on government decisions, changing testing policies as well as events with impact on epidemiological dynamics such as holidays and sudden outbreaks(such as thin peak of new infections in June affected mostly workers of the meat industry). We also assumed additional time points with changes determined by BIC.
Again, for the sake of parsimony, we assume that the relative infection intensities of asymptomatic ( ( )) and symptomatic subjects ( ( )) are the same, hence, respective proportionality , is constant and estimated during model fitting. Thus, ( ) is defined as follows: 2. Likewise, rates toward critical disease states and deaths are also assumed to vary through the course of the epidemic due to changes in testing policy resulting in different percentages of unreported cases and asymptomatic subjects, changes in age distribution of infected subjects, improvement of patient care due to new treatment options and due to possible over-stretched medical resources (not the case in Germany but other countries). In our model, this is also accounted for phenomenologically by assuming the probabilities and as time-dependent input parameters. We assume that both functions are step functions: where , are empirical dates and The partitions of pcrit and pdeath are assumed independent. Respective numbers of jumps and can differ.
The time point t = 0 corresponds to the 3 March 2020.
2. Likewise, rates toward critical disease states and deaths are also assumed to vary through the course of the epidemic due to changes in testing policy resulting in different percentages of unreported cases and asymptomatic subjects, changes in age distribution of infected subjects, improvement of patient care due to new treatment options and due to possible over-stretched medical resources (not the case in Germany but other countries). In our model, this is also accounted for phenomenologically by assuming the probabilities p crit and p death as time-dependent input parameters.
We assume that both functions are step functions: where T pcrit,j N crit j=1 are empirical dates and α crit,j N crit j=1 the respective relative changes of p crit . Both, T pcrit,j N crit j=1 as well as α crit,j N crit j=1 are parameters to be estimated. The initial value of p crit is p crit,0 . Functions χ [t j ,t j+1 ) (t) are indicator functions being 1 in the interval T pcrit,j , T pcrit,j+1 and 0 else. The step functions for p death (t) and p death,S (t) are defined analogously: The partitions of p crit and p death are assumed independent. Respective numbers of jumps N crit and N death can differ.
In order to find an optimal tradeoff between parsimony and goodness of fit we calculated Bayesian Information criteria (BIC) for different partition numbers N tr , N crit , and N death and chose partitions minimizing BIC.
When new data become available, we attempt to update the numbers of partitions every two weeks by considering adding a new break-point within the last month. The time point as well as the corresponding jump value are considered as two new parameters. We added a new break point only if it improves BIC after the updated parameters estimation.

Appendix C. Output Layer
We here describe, how the state parameters of the hidden SECIR model are linked with data via the output layer of the IO-NLDS.
Modeling of daily registered infected cases I S,M : The total number of daily registered infected cases I S,M is coupled to the efflux of the first asymptomatic compartments I A,1 and I A, 1 Mu toward symptomatic compartments multiplied with P S,M . However, we assume a delay of the registered vs. reported cases by introducing an empirical distribution of the reporting delay: When a person has a positive PCR SARS-CoV-2 test result at a date d 1 , this will result in a registered case at a later date d 2 . The difference d = d 2 − d 1 is the reporting delay and is assumed log-normally distributed. This distribution of delays is determined on the basis of data provided by the Robert-Koch-Institute from the period 27 April 2020 to 13 November 2020. The parameters of this distribution are derived by minimizing the Kullback-Leibler divergence between the parametric representation and the empirical distribution. Results are displayed in Figure A1. Modeling delay of death reporting: In contrast to the newly infected cases, neither information of delays in COVID-19 associated death reporting nor actual dates of deaths were available to us. Therefore, we used a data model proposed by Delagdo et al. [45]. In detail, we assumed that the delay is normally distributed with an average of 7.14 days and a standard deviation of 4 days.
Since we consider time as an integer, we discretize this normal distribution by the approximation DRD(d) = N 7.14,4 (d) ∑ 100 i=1 N 7.14,4 (i) for integers d ≤ 100 and 0 else, where N is the Gaussian distribution function with mean 7.14 days and standard deviation of 4 days, i.e., we neglect delays larger than 100 days. Using this approximation, we derive the actual number of new deaths at time point t: Here, D r (t) is the number of reported new deaths at time t. The function D a (t) is linked to our compartment D.

Appendix D. Parameter Estimation
Free parameters of the model are determined by minimizing the negative log-likelihood function of observed data. The negative likelihood is constructed in analogy to [47]. It constitutes of the sum of three components: The terms nLL pri i and nLL resid i correspond to prior constraints of parameters and to the residual errors of the data as explained below in detail. The term Constr is a penalty term to keep values in eligible ranges or orders (see "Penalization"). We assume independence between parameters throughout.
Parameter distributions and transformations: Most of the parameters are confined to certain ranges. During estimation (with possible prior constraints), we transform these parameters to the space of real numbers. We assume that these transformed values are normally distributed during Markov-Chain Monte Carlo (MCMC) sampling (see below). To ensure this, parameters confined to a finite interval (a,b) are transformed by the logitfunction. Parameters with positive values are transformed by a log-normal transformation. Thus, where ϕ s is the s-th parameter and ψ s is the respective transformed parameter and N par is the total number of parameters to be estimated. The negative likelihood contribution of the priors nll pri i is defined as follows: where δ s equals 1, if a prior is assumed for the s-th parameter and 0 otherwise. The prior information is represented by the "best value" ψ pri s and an uncertainty expressed as standard deviation of possible values ω pri,s . We assume that parameter estimates are random variables normally distributed around their respective prior values. Thus, Prior "best values" and ranges of parameters are provided at Table 2 of the main paper. The uncertainties ω pri,k are set to 2 for all parameters. This heuristic setting is based on a tradeoff between avoidance of overfitting including implausible parameter values and good data fitting properties.
Penalization: We penalize with a high value of 10 8 in cases when times of nonpharmaceutical interventions are either too close (closer than 3 days) or non-monotonic. In the same way, we penalize too high dynamical p death values (more than 0.66) by multiplication of max(p death -0.66,0) with 100.
Residual errors of observed vs. predicted data: We fit data for daily registered cases, cumulative registered cases, deaths, cumulative deaths, and ICU occupation as explained in sub-section "output layer" and the methods section of the main paper. The respective term of the negative log-likelihood nll resid i corresponds to the residual errors of these data. Thus, maceutical interventions are either too close (closer than 3 days) or non-monotonic. In the same way, we penalize too high dynamical pdeath values (more than 0.66) by multiplication of max(pdeath-0.66,0) with 100.
Residual errors of observed vs. predicted data: We fit data for daily registered cases, cumulative registered cases, deaths, cumulative deaths, and ICU occupation as explained in sub-section "output layer" and the methods section of the main paper. The respective term of the negative log-likelihood corresponds to the residual errors of these data. Thus, where represents the output layers ("dY" corresponds to daily counts, while "Y" corresponds to cumulative counts). Subscript S denotes simulation results, while D corresponds to the data. We sum the negative log-likelihoods of the three outputs considered (infected subjects, critical cases and deaths). Thus, represents one of these three entities x with number of data points at time points , (j = 1,…,Nx) and residual errors ax. We introduce weights for the cumulative terms as compared with the daily counts and set it to 0.2. The cumulative terms were introduced to avoid biases of cumulative data occurring after fitting daily data only. Cumulative data for ICU occupation were not fitted, i.e., = 0. The parameter corresponds to the power transformation used to compare model and data. In the present model version, it is set to 0.5. It constitutes a tradeoff between fitting precision of large and small numbers. All weights and were set to 1. Thus, we assume that for each output and for each data point the entities The algorithm to minimize the negative log-likelihood is explained in the next section. Differences of estimated values and their respective priors can be tested by calculating Z-scores , .

Appendix E. Algorithm for Parameter Estimations and Prediction Sampling
Due to nonlinearity of (Equations (A8), (A11) and (A12)), parameters and residual errors θ cannot be estimated simultaneously. For such situations an expectation-maximization (EM) algorithm was proposed by Dempster et al. [58]. This algorithm is a widely applied approach for the iterative computation of maximum likelihood (correspondingly minimum of negative log-likelihood) estimates in incomplete-data statistical problems. In detail, the random parameters = are considered as non-observed data, while observed data y in our case are defined as follows: where Y out represents the output layers ("dY" corresponds to daily counts, while "Y" corresponds to cumulative counts). Subscript S denotes simulation results, while D corresponds to the data. We sum the negative log-likelihoods of the three outputs considered (infected subjects, critical cases and deaths). Thus, Y out represents one of these three entities x with number of data points N x at time points t j,x (j = 1, . . . ,N x ) and residual errors a x . We introduce weights we cumul for the cumulative terms as compared with the daily counts and set it to 0.2. The cumulative terms were introduced to avoid biases of cumulative data occurring after fitting daily data only. Cumulative data for ICU occupation were not fitted, i.e., we ICU = 0. The parameter tr Y out corresponds to the power transformation used to compare model and data. In the present model version, it is set to 0.5. It constitutes a tradeoff between fitting precision of large and small numbers. All weights we dY out and we Y out were set to 1.
Thus, we assume that for each output and for each data point the entities dY

Appendix E. Algorithm for Parameter Estimations and Prediction Sampling
Due to nonlinearity of (Equations (A8), (A11) and (A12)), parameters ψ and residual errors θ cannot be estimated simultaneously. For such situations an expectation-maximization (EM) algorithm was proposed by Dempster et al. [58]. This algorithm is a widely applied approach for the iterative computation of maximum likelihood (correspondingly minimum of negative log-likelihood) estimates in incomplete-data statistical problems. In detail, the random parameters ψ = {ψ s } N par s=1 are considered as non-observed data, while observed data y in our case are defined as follows: Complete data of the model is (y, ψ). The unknown residual errors θ describe the uncertainty of parameters ψ.
Therefore nLL(y, ψ; θ) is a marginal negative log-likelihood likelihood. The complete likelihood nLL is defined as follows: The EM algorithm minimizes nLL(y; θ) iteratively: At the k-th iteration of EM, the expectation step computes the conditional expectation of the complete negative log-likelihood Q k (θ) = E(nLL(y, ψ; θ)|y, θ k−1 ) by generating ψ (k) based on previous estimates θ k−1 , and the maximization step computes the value θ k maximizing Q k (θ). The EM sequence (θ k ) converges to a stationary point under general regularity conditions [58].
In nonlinear cases, the expectation step cannot be performed in a closed form. Therefore, we applied the Stochastic Approximation algorithm of EM (SAEM). SAEM is a maximum likelihood estimator of the population parameters [48] based on stochastic integration of marginal probabilities without likelihood approximation such as linearization or quadrature approximation or sigma-point filtering [17]. Our implementation is inspired by and is very similar to that of earlier versions of Monolix (Lixoft) software (http://lixoft.com/, accessed on 16 November 2018) The stochastic approximation version of standard EM algorithm (SAEM) proposed by [48] replaces the usual E-step at an iteration k by a stochastic procedure as follows: 1.
Simulation step: draw m k realizations of ψ (k) = ψ (k) s N par s=1 from the conditional distribution p(·|y; θ k ) using MCMC algorithm.

2.
Stochastic approximation: update Q k (θ) where γ k is a decreasing sequence of positive numbers.

1.
Our stochastic approximation step is an improved version of the stochastic approximation of the integration of marginal distribution on the multidimensional domain Ω of possible parameter values: In analogy to Monolix software, we selected γ k as follows: We choose K 1 equal to 4 and run the algorithm until convergence with a tolerance 0.1% of estimates of population parameters (see below).

3.
We performed MCMC sampling 4000 times at each stage with a burn-in phase of 1000 steps. Thus, m k = 3000.
Exact estimates of different components of θ k are: Therefore, respective stochastic approximations and maximizations θ k are as follows: In the same way, the respective terms for the cumulative data approximations are derived.

Appendix F. MCMC Algorithm for the Expectation Step
Markov chain Monte Carlo (MCMC) methods comprise a class of algorithms for sampling from a probability distribution [59]. By constructing a Markov chain that has the desired distribution as its equilibrium distribution, one can obtain a sample of the desired distribution by recording states from the chain. It is well-known that a proper choice of a proposal distribution for MCMC methods is a crucial factor for convergence of the algorithm [60]. For the sake of increasing the acceptance rate, a number of adaptive Metropolis (AM) algorithms were proposed by different groups. Here the proposal distribution is learned along the process using the full information cumulated so far. We implemented the adaptive MCMC version with Gaussian proposal distribution described in [60] as well as adaptive incremental Mixture MCMC [56] called AIMM, which we modified slightly. Strictly speaking, these methods are not really Markov chains, because proposal distribution of the next step depends on all preceding states → X t t 0 rather than only the previous one. The algorithm of Haario et al. is simpler and it assumes the existence of a global minimum of nLL. In contrast, the algorithm of Maire et al. could be useful for cases when the nLL has a complex topology due to overfitting.
Let π denote the target distribution (i.e., negative log-likelihood) given by Equation (A8). At each iteration step the new parameter vector Y is generated by a transition kernel representing the proposal distribution. This candidate vector is accepted with probability The transition kernel of Haario's MCMC version is an empirical covariance matrix of previous samples stabilized by an identity matrix multiplied by a small number ε: . Here, we choose ε = 0.0001. This parameter is required to ensure ergodic property of the Markov chain. When k > 1, we added samples from the previous step to the covariance matrix: At each iteration k we used a sample from the previous iteration providing a small value of nLL as starting point.
In the AIMM, the proposal distributions Q t are mixtures of multivariate normal distributions. Roughly spoken, this is a generalization of Haario's algorithm when multiple local minima of the nLL exist in few clusters. The candidate vector is accepted with probability This algorithm minimizes discrepancies between proposal and target probability i.e., a sequence {Q t } converges to π by approximating it through mixtures of multivariate normal distributions. The elements of this series Q t are defined as follows: where M t is the number of components at the iteration t. The elements {ϕ 1 , · · · , ϕ M t } represent the incremental mixture components, {β 1 , · · · , β M t } are the respective weights. After the choice of the r-th component, a random parameter vector Y is generated around the r-th mean according to the r-th covariance matrix as in Haario's algorithm. If Y is accepted, it becomes → X t . → X t can either stay in the r-th cluster or give origin for the new cluster ϕ M t with M t = 1 + M t−1 . A new cluster is created when the match of ϕ M t to the r-th cluster is insufficient based on Mahalanobis distance. If → X t stays in the r-th cluster, it updates the r-th covariance matrix in a similar way as in Haario's algorithm [60], Equation (A26).
We here modified the conditions for new cluster formation compared to [56] as follows. In our algorithm, a new cluster is formed when one of the following conditions hold:

•
The Mahalanobis distance of → X t to the cluster from which it was generated is less than 0.025 or larger than 0.975. That is → X t diverges significantly from the current multivariate normal distribution of the r-th cluster • π → X t is significantly larger than π of the current cluster center. That is → X t does not correspond to the local maximum of π in the neighbourhood of the r-th cluster.
If one of the above conditions holds, → X t becomes the center of a new cluster. The respective Gaussian component is the covariance matrix of the r-th (i.e., previous) cluster. This matrix will be further updated every time when new members of the new cluster are accepted in future proposals.
The weights {β 1 , · · · , β M t } are proportional to π of the respective cluster centers to the power of γ, where γ is a positive number less than 1. All weights are updated every time when a new cluster emerges.
In summary, AIMM accepts proposals with discrepancies to the target distribution. As a consequence, proposal distributions are multivariate normal mixtures. Every cluster's mean is a local maximum of π. Sampling of proposal distributions from clusters depends on π. New clusters emerge when an accepted proposal either significantly diverges from the respective cluster's probability or when a significantly better optimal value is found in this cluster.
After thorough comparison of adaptive MCMC and the adaptive incremental mixture MCMC, we found the latter to be superior. Higher values of π were found in a shorter time. It also generates higher acceptance rates (0.2-0.3 versus 0.1) and finds more alternative solutions. We therefore used this method for our parameter estimations.
We applied Geweke convergence diagnostics for Markov chains [61] as implemented in the R-package coda (https://cran.r-project.org/web/packages/coda/coda.pdf, accessed on 01 October 2020). This method is based on a test for equality of the means of the first and last part of a Markov chain (by default the first 10% and the last 50%). If the samples are drawn from a stationary distribution of the chain, the two means are equal and Geweke's statistic has an asymptotical standard normal distribution. The test statistic is a standard Zscore: the difference between the two sample means divided by its estimated standard error. The standard error is estimated from the spectral density at zero taking autocorrelation into account. The Z-score is calculated under the assumption that the two parts of the chain are asymptotically independent. We applied this diagnostic for the nLL resulting from the last run of the implemented adaptive MCMC algorithm, resulting in a chain of 2620 steps. We considered default fractions of the chain, i.e., 0.1 and 0.5 of the beginning and from end of chain, respectively. The resulting Z-score is −0.2536, corresponding to p-value of 0.4, i.e., no deviations of means were detected suggesting that a stationary distribution is achieved. Figure A2 (generated by function geweke.plot from the coda package) shows the development of Geweke's Z-score when successively larger numbers of iterations are discarded from the beginning of the chain. The Z-score remains always in the 95% confidence interval, suggesting a successful convergence.

Appendix G. MCMC Simulation for Prediction and Controlling Goodness of Fit
The estimates of residual errors are determined at the last step and are used for MCMC sampling of parameters ψ. The resulting means and standard deviations are considered as respective average estimates and their standard errors. Simulations of these parameter samples provides a set of alternative predictions. From these, we collected the best fitting solution, the average solution and confidence intervals for different confidence limits α.

Appendix H. Justification of Prior Parameters and Ranges
We here provide justifications of assumed prior values and parameter ranges. Details of parameters definition and fitting can be found in Appendices B and D, respectively.
Initial influx of people per day Influx: The initial influx was estimated from the data without prior assumptions to a value of 6937 people per day in order to initialize the simulation. Later, the parameter is no longer relevant for simulation outcomes.
Infection rate through asymptomatic subjects per day r 1 : This infection rate was estimated from the data without prior assumptions. It represents the basic transmission probability of the SARS-CoV-2 virus from an asymptomatic infectious person to a susceptible contact.
Infection rate through symptomatic subjects per day r 2 : This infection rate was estimated from the data without prior assumptions. It represents the basic transmission probability of the SARS-CoV-2 virus from a symptomatic infectious person to a susceptible contact.
Relative infection intensity of asymptomatic subjects per day b 1 (t): The infectivity of asymptomatic infected subjects was assumed as a time-dependent step-function due to changing NPIs/contact behavior and other factors influencing infection probabilities. Steps were estimated from the data without prior assumptions.
Relative infection intensity of symptomatic subjects per day b 2 (t): The infectivity of symptomatic infected subjects was assumed as a time-dependent step-function due to changing NPIs/contact behavior and other factors influencing infection probabilities. Steps were estimated from the data without prior assumptions.
Ratio of b 1 (t) and b 2 (t) (rb 1,2 ). We assumed a fixed ratio of the infectivities of asymptomatic and symptomatic infected subjects for the sake of parsimony. The ratio was estimated from the data without prior assumptions.
Fixed dates for updates of infectivity functions: We used several fixed dates of changes in infectivity functions due to known changes in NPIs, testing policy or outbreaks. Note that even fixed time points were checked for necessity to assume changes in infectivity for the sake of parsimony, i.e., respective steps were only assumed if significantly improving model fit. The first three fixed time-points, tr 1 (10 March 2020), tr 2 (15 March 2020), and tr 3 (22 March 2020) reflect German governmental interventions including regulation of the size of public events, travel restrictions, and contact restriction. Fixed time-points tr 6 (30 April 2020), tr 7 (7 May 2020), and tr 8 (21 May 2020) reflect German governmental interventions related to the step-wise relaxation of NPIs, in particular regarding leisure sports, contacts, and schools. Time point tr 17 (2 November 2020) reflects governmental NPIs in response to the German second wave, including restrictions of public life and social contacts, also referred as "soft lockdown". Time point tr 21 (16 December 2020) reflects further stricter governmental NPIs in response to the ongoing increase of the German second wave, also referred as "hard lockdown", strongly limiting public and private contacts including school closures. Finally, time point tr 28 (23 February 2021) reflects release of many governmental NPIs in response to the decline of the German second wave.
Transit rate for compartment E (latent time) r 3 : The transition rate r 3 for the compartment of exposed subjects is the inverse of the latent time, i.e., the time being infected but not yet infectious. The mean of the prior distribution for the latent time was set to 3 days and the minimum and maximum of the distribution was set to 2 and 4 days, respectively, in accordance with previous reports [10]. Note that minimum and maximum of a parameter's distribution in this section always refer to the distribution of the mean of the parameter, not of the distribution of the parameter itself. Further justification of this parameter is discussed in the following when considering the rate r 4b .
Transit rate for asymptomatic sub-compartments r 4 : The transition rate r 4 for the asymptomatic infectious compartment to the recovered compartment is a third of the inverse of the time being asymptomatic and infectious, as this compartment is split into three sub-compartments. The mean of the prior distribution of r 4 was set to 3/5 per day and the minimum and maximum of the distribution was set to 3/10 and 3/4 per day, respectively. These values are based on general considerations regarding timelines of the germinal center reaction [25] and further supported by reports from the literature estimating relevant infectiousness periods in general or asymptomatic/mildly symptomatic COVID-19 patients as in between 3.5 and 9.5 days [22][23][24].
Rate of development of symptoms after infection r 4b : The inverse of this rate is equal to the time from being infectious to start of developing symptoms. The mean of the prior distribution of r 4b was set to 1 2 .5 per day and the minimum and maximum of the distribution was set to 1/5 and 1/1 per day, respectively. This is in line with previous reports [10,26,27]. Note that the serial interval, i.e., the average time between successive cases in a chain of transmission is composed of two parameters of our model. In detail, the serial interval is the sum of the average latent time (1/r 3 ) and half of the average time being infectious when assuming random occurrence of subsequent infections during time of infectiousness. Exemplarily, if the serial interval would be considered in a scenario where symptomatic individuals are immediately and effectively quarantined, the serial interval would be 1/r 3 + 0.5*1/r 4b . The serial interval was estimated by the RKI [28] to have a median of 4 days (interquartile range 3-5 days) based on the literature [18][19][20][21], which is in accordance with our choices for r 3 and r 4b . However, the serial interval (and other parameters like the time being infectious) are to some extent also time dependent, reflecting e.g., behavioral changes. Although we do not model a time dependence for these specific parameters, our model can, to a certain extent, cope for this by data-driven adaptation of other time-dependent parameters such as b 1 and b 2 .
Probability of developing symptomatic disease after infection p symp : This probability was estimated from the literature, reporting a percentage of symptomatic COVID-19 cases in between 55% and 85% [37][38][39]. We used a percentage of 50% as mean of the prior distribution, accounting for the fact that minor symptoms are frequently ignored or considered as symptoms of a common cold. Minimum and maximum was set to 0.3 and 0.8, respectively.
Transit rate of symptomatic sub-compartments r 5 : The transition rate r 5 for the three symptomatic sub-compartments towards recovery is a third of the inverse of time being symptomatic and infectious. The mean of the prior distribution of r 5 was set to 3/2.5 per day and the minimum and maximum of the distribution was set to 3/7.5 and 3/1.5 per day, respectively. These values are based on the assumption that symptomatic and asymptomatic subjects are similar with respect to time of contagiousness. Hence, values of the distribution of r 5 equal that of r 4 subtracted by the mean value of r 4b .
Rate of development of critical state after becoming symptomatic r 6 : The inverse of this rate is assumed equal to the time of developing a critical state after being infectious and symptomatic. The mean of the prior distribution of r 6 was set to 1/5 per day and the minimum and maximum of the distribution was set to 1/7 and 1/4 per day, respectively, according to previous reports [10,[29][30][31]. Note that the probability of people becoming critical is affected by the function p crit .
Probability of becoming critical after developing symptoms p crit : This probability is assumed as a time-dependent step-function reflecting for example changing age-distributions of infected subjects or treatment efficacy. Steps were estimated from the data within the range of 0 to 1 without assuming a specific prior. The initial value of p crit,0 was estimated as 0.075, which is within the range of reported values [10,62].
Transit rate for critical state sub-compartments r 7 : The transition rate r 7 for the critical state sub-compartments is a third of the inverse of the time treated in intensive care unit (ICU) for survivors, as the critical compartment is also split into three sub-compartments. The mean of the prior distribution of r 7 was set to 3/17 per day and the minimum and maximum of the distribution was set to 3/35 and 3/8 per day, respectively. These values are informed by previous reports focusing on data of 35 to 79-year-old patients, the most frequent population in ICU [10,[32][33][34].
Death rate of patients in critical sub-compartment r 8 : This transition rate represents the rate from the first ICU sub-compartment to the death compartment. It is the inverse of the time of patients in ICU that passed away. The mean of the prior distribution of r 8 was set to 1/8 per day and the minimum and maximum of the distribution was set to 1/14 and 1/6.5 per day, respectively. This reflects the shorter time in ICU for patients with fatal disease outcome informed by previous reports [29,35,36]. The number of people with fatal disease course is affected by two additional parameters p death and p death,S explained below.
Probability of death after becoming critical p death : This is the probability of death for patients at ICU. It is assumed as a time-dependent step-function estimated from data. Values are restricted within the range 0 to 1 without specific prior assumptions. Changes in time reflect for example changes in age-composition of ICU patients as well as changes in treatment regimens. The initial value is p death,0 = 0.118.
Probability of death after developing symptoms without becoming critical p death,S : To reflect COVID-19 related deaths outside of ICU (especially relevant for the oldest age-groups [32]), we introduced the probability p death,S of transitioning from the second symptomatic sub-compartment to the death compartment. This probability was estimated from the data p death,S = 0.0448.
Fraction of unreported cases p S,M : For the fraction of infected cases that are symptomatic but not reported, we used a prior distribution with a mean of 0.5, a minimum of 0.3 and a maximum of 0.9. This choice was informed by studies of SARS-CoV-2 seroprevalence in Germany [40,41]. Note that the total percentage of unreported infected people is 1-p S,M ·p symp according to the definition of p symp .
The factor mur is multiplied to r 1 and r 2 reflecting higher infectivity of the B.1.1.7 variant compared to the previous variants. This parameter was estimated from sequencing data reporting the dynamics of the increase of variant B.1.1.7 in the UK, Denmark, Belgium, Suisse, and the United States, available from https://github.com/tomwenseleers/newcovid_ belgium/, accessed on 13 April 2022) and Germany, available from "Mutationstracking-Projekt von Sven Schmidt" at https://tinyurl.com/36xnmxat, accessed on 17 May 2022). Thereby, mur was calibrated to match the observed average dynamic of the increase of B.1.1.7 across countries resulting in a value of mur = 1.7. Table A1. Time points of changes in infectivity and respective steps. We used fixed (known due to Governmental decisions or random events) and estimated time points of NPI/contact behavior changes and events and respective changes in infectivity of asymptomatic subjects. We provide estimates and relative standard errors of the infectivity. For estimated time points, we also provide the respective standard error (last column).

Number
Type  Table A3.
Step functions of p crit and p death . We present estimates for the single steps of the functions p crit and p death at the specified dates and respective standard errors. We also provide the standard error of the estimated time point (last column).   Table A5. Parameter estimates and comparison with average priors. We present estimated parameters of the SECIR model and initial conditions of control parameters and their respective standard errors for Germany. We also perform a formal comparison of estimates and expected priors using t-test.  Table A6. Time points of changes in infectivity and respective values for Saxony. We used fixed (known due to Governmental decisions or random events) and estimated time points of changes of NPI/contact behavior and events and respective changes in infectivity of asymptomatic subjects. We provide estimates and relative standard errors of the infectivity starting with the date mentioned (3 to 5 column).  Table A7.

Numbers
Step functions of p crit and p death for Saxony. We present estimates for the steps of the functions p crit and p death at the specified dates and respective standard errors for Saxony. We also provide the standard error of the estimated time point (last column).   Table A8. Parameter estimates and comparison with average priors for the parameter settings for Saxony. We present estimated parameters of the SECIR model and initial conditions of control parameters and their respective standard errors for the parametrization of the epidemic in Saxony. We also perform a formal comparison of estimates and expected priors using t-test.