Early stage COVID-19 disease dynamics in Germany: models and parameter identification

Since the end of 2019 an outbreak of a new strain of coronavirus, called SARS-CoV-2, is reported from China and later other parts of the world. Since January 21, World Health Organization (WHO) reports daily data on confirmed cases and deaths from both China and other countries (www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports). The Johns Hopkins University (github.com/CSSEGISandData/COVID-19/blob/master/csse_COVID_19_data/csse_COVID_19_time_series/time_series_COVID19_confirmed_global.csv) collects those data from various sources worldwide on a daily basis. For Germany, the Robert-Koch-Institute (RKI) also issues daily reports on the current number of infections and infection related fatal cases (www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Situationsberichte/Gesamt.html). However, due to delays in the data collection, the data from RKI always lags behind those reported by Johns Hopkins. In this work we present an extended SEIRD-model to describe the disease dynamics in Germany. The parameter values are identified by matching the model output to the officially reported cases. An additional parameter to capture the influence of unidentified cases is also included in the model.


Introduction
In December 2019, first cases of a novel pneumonia of unknown cause were reported from Wuhan, the seventh-largest city in China. In the meantime, these cases have been identified as infections with a novel strain of coronavirus, called SARS-CoV-2 and the disease it causes is called coronavirus disease 2019 . At the beginning of January 2020, the virus spread over mainland China and reached other provinces. Increased travel activities due to the Chinese new year festivities supported the expansion of the infection. Since 21 January, WHO's daily situation reports contain the latest figures on confirmed cases and deaths, see [1].
The first COVID-19 case in Germany was reported in late January 2020 in a company close to Munich, Bavaria. Later cases were imported by travelers from China, Iran or Italy as well as tourists returning from ski holidays in the Austria and Italy. By 1 March 2020 more than 100 cases were reported in Germany and since than the number of cases began to rise exponentially, see Fig. 1. The first deaths were reported on 9 March 2020 [2,3].
By 16 March 2020 the federal government introduced first measures to reduce the spread of the disease: Schools, kindergartens and universities were closed. On 22 March these measures were tightened by implementing a national curfew and contact ban. People are advised to stay at home, leaving only for work related activities, necessary shopping, medical treatment or sports. All this should not be done in groups of more than two persons if they do not belong to the same household [4].
Our work is based on the data reported by Johns Hopkins University [5]. We refrain from using the official data from the Robert-Koch-Institute [2], since they suffer from a delay by several days due to the more complicate way of aggregating those data. For a detailed explanation of the difference between the data reported by Johns Hopkins and the Robert-Koch Institute we refer to the information given on the webpage of the Robert-Koch-Institute, see [6]. Johns Hopkins University continuously collects the data from internet queries at various sources (local health authorities, newspapers, etc.) whereas the Robert-Koch-Institutes collects the data that are reported for the local health authorities to the district level, then state level and finally aggregates them to the federal statistics. Hence these data lag several days behind the ones collected by Johns Hopkins University.
The paper is organized as follows: In Sect. 2 we describe the model and the parameter identification problem. Our models consists of three variants of a five compartment SEIRD-system without demographic terms, where the transmission rate is either fixed (1a)-(1e) or time-dependent (3a)-(3e) and (4a)-(4e). The fatalities are either described by an ODE, see models (1a)-(1e) and (3a)-(3e), or via a delay term in model (4a)-

Figure 2
Transmission diagram for the basic SEIRD-model (1a)-(1e). The artificial compartment C contains all infected cases, i.e. current active infections, recovered and deaths (4e). In the parameter estimation problem, we determine the transmission rate, detection rate and lethality together with the initial values for the exposed and infected compartment. In Sect. 3 we discuss the sensitivity of our model with respect to detection rate. Section 4 is devoted to the adjoint equations used for solving the optimization problem. The simulation results are presented in Sect. 5. Here we do compare the results obtained from the three models presented in Sect. 2.

Mathematical model
To model the dynamics of the spread of COVID-19 incidences, we propose a hierarchy of SEIRD models. For details regarding the original SIR-and SEIR-model we refer to classical works on mathematical epidemiology, e.g [7]. For our basic SEIRD-model, the total population of Germany with N ∼ 83.000.000 individuals is subdivided in to susceptibles S, exposed E, infected I, recovered R and deaths D. The susceptibles constitute the reservoir of persons that are not yet infected with SARS-CoV-2. After infection susceptible become exposed meaning that they already carry the virus but are not yet infectious. With a rate ϑ exposed individuals become infectious and transmit the virus with rate β to susceptibles. An infected individual loses infectivity with γ and has a probability μ of dying due to the disease [8]. Figure 2 shows the transmission structure. By C we denote all infected cases, independent of their current status. This artificial compartment is later on used to compare with the total number of registered cases reported by Johns Hopkins or RKI.
The resulting system of ordinary differential equations (ODE) for the above described SEIRD-model reads as The starting time t 0 is chosen as 1 March and the initial conditions for the recovered and dead compartment are assumed to be zero, since in Germany the first COVID-19 In the sequel, we will also consider two refined versions of the above basic model. At the onset of the disease, the numbers of exposed, infected, recovered and dead are still small and the number of susceptibles is approximately equal to the entire population N . In this setting, the EI-part of the model reduces to The maximal eigenvalue λ of this linear system determines the initial growth rate and is given by and the doubling time T 2 equals To include this into the basic model (1a)-(1e), we also consider an alternative model for the transmission rate β: We assume β as a piecewise constant function on the time intervals prior to any measures, (until 15 March), after school closings (between 16 and 22 March) and after the contact ban (after March 22) The resulting time-dependent SEIRD-model reads as Setting β := β 0 = β 1 = β 2 , the time-dependent model reduces to the basic one. In order to validate our models and to identify the parameters involved therein, both the registered number of infections and the registered number of COVID-19 related deaths are important indications. The number of registered deaths is probably considerably more reliable, since the number of registered infections depends on the number of tests conducted and the dark figure of undetected, mostly asymptomatic cases, is assumed to be remarkably large [9]. We will discuss this point later in more detail. In the previous basic or time-dependent SEIRD-model, the actual increase of the disease related deaths dD dt is assumed to be proportional to the current number of infected persons. The Robert-Koch-Institute specifies an average of 10 days between the onset of symptoms and admission to the intensive care unit [10]. Therefore, we assume τ = 14 for the time between the onset of infectiousness and death. In order to include this time lag into our model, we introduce a delay-term into the time-dependent model and obtain the final delayed time-dependent model: Note, that for solving this delay differential equation (DDE) we need an initial history of the infected compartment, i.e. values I 0 (s) for t 0τ ≤ s ≤ t 0 .
In all the three models, the parameters ϑ = 1/2 [days -1 ], γ = 1/10 [days -1 ] are assume to be fixed and resemble a latency period of 2 days and a recovery period of 10 days, see [2, Situation report 31 March 2020].
The parameters in the transmission rate, i.e. β, or β 0 , β 1 , β 2 the lethality μ and the initial values E 0 , I 0 resp. the initial history I 0 (s) for the exposed and infected compartment are yet unknown to us. We will identify them together with the detection rate δ by matching the model output to the given data. The detection rate δ corresponds to the fraction of infected individuals which are positively tested for SARS-CoV-2 and hence appear in the official recordings. Various sources speculate that this detection rate is in the order of magnitude of 10-20% meaning that the true number of infected 5-10 times larger than the number published in the official statistics, see [9].
To match the model output and the reported data we use a least-squares approach. Let u = (β, δ, μ, E 0 , I 0 ) resp. u = (β 0 , β 1 , β 2 , δ, μ, E 0 , I 0 ) denote the unknown model parameters to be determined. Furthermore, let Y (t) and Z(t) denote the data for the cumulated infected and dead cases at time t reported by Johns Hopkins University. The deviation between the model and the data is measured by the cost functional where f 2 1.2 · 10 11 and Z 2 L 2 6.5 · 10 8 , hence ω 1 c 1 · 185. The cumulated infected Y , i.e. total positive tests, are to be matched in the SEIRD-model to those individuals who had been infected until time t, i.e. the sum of the infected I, recovered R and deaths D. To account for the uncertainty in the true number of infected and recovered cases, we multiply both compartments by the detection rate δ, which is itself part of the parameters to be identified. For the deaths we assume no undetected cases. By T Fit we denote the time horizon used for the comparison between the model and the data. The regularization term ω 2 u 2 is included to ensure the convexity of the cost-functional. The weighting parameters c 1 , c 2 and hence ω 1 , ω 2 > 0 allow to balance the contributions from the least squares error in the fatalities and from the size of the parameter values themselves to the least squares error in the infected cases. The weight c 1 for the fatal cases allows to compensate the different order of magnitude between the infected cases and the fatal cases, typically c 1 2-3 leading to ω 1 500. The weight c 2 is chosen small, such that the overall cost functional is still dominated by the least square fit between the model output and the given data.
The parameters u * themselves are obtained from minimization problem

A few analytical considerations
Due to the absence of demographic terms, our basic model (1a)-(1e) does not allow other equilibria besides the trivial disease free equilibrium X 0 = (N, 0, 0, 0, 0). Since we focus only on the short-time behavior of the epidemics, demographic terms are excluded and equilibria do not play any important role. An important issue is the question of wether we can identify the detection rate and lethality during the take-off period of the epidemics? The only data available for parameter identification are the total number of registered cases C = I + R + D and the deaths D. The total registered cases heavily depend on the number of tests conducted. If a person is infected, but not tested, this person will not appear in the official statistics. Hence, there is a presumably large dark figure in the officially recorded data. Our model parameter δ takes this into account. The other, maybe more reliable, available data are the recorded deaths.
Here we may assume that all COVID-19 related deaths are diagnosed and hence there is no dark figure in the D-compartment. A recent analysis by the Federal Statistical Office on the excess mortality in Germany for March and April 2020 confirms this assumption, see [11]. For other countries this assumption might be questionable, since they suffered from major COVID-19 outbreaks in care homes that did not enter the official statistics, e.g. in the UK, see [12].
However, one scenario could be possible. A large dark figure in the entire cases, i.e. a small detection rate δ and a very small lethality could result in the same or at least similar observed data as a moderate or even small dark figure and hence large detection rate δ combined with a higher lethality rate. In that setting a simultaneous identification of both, the detection rate δ and the lethality μ could be difficult due to their counteracting effects.
In order to investigate this scenario, we consider the simultaneous effect of the detection rate δ scaling both the initial values of the E and I compartment to account for undetected cases together with a lethality δμ. Removing the S-compartment by setting S = N -E -I -R -D, the basic SEIRD-system (1a)-(1e) reads as The sensitivities Σ E := ∂ δ E and Σ I , Σ R , Σ D of the solution with respect to the detection rate satisfy the system , Σ I = ϑΣ Eγ Σ I , Σ I (t 0 ) = -I 0 /δ 2 , ( 7 b ) To illustrate these findings, we consider a linearization of a simplified SIR-model during the initial phase of the epidemics. We neglect the exposed compartment and assume that at the initial phase, the number of susceptibles is approximately equal to the entire population. Hence we get the linear system In this linearized setting, the approximationD for the dead compartment is independent of the detection rate δ.
From the graphs in Fig. 4 one can conclude, the a significant dependence of the detected or dead compartment C resp. D is given only after the initial take-off period of the epidemic. In the setting of Germany, this implies, that during the month of March a reliable identification to the detection rate might not be possible.

Adjoint equations and optimization
In order to solve the minimization problem (6a)-(6b), we use the adjoint equations, for details see [13,14]. We introduce the Lagrangian Here z = (z S , z E , z I , z R , z D ) denotes the adjoint functions to the state variable x = (S, E, I, R, D) and g(t, x, u) denotes the right hand side of the ODE resp. DDE system. The gradient of L with respect to the unknown parameters u is given by , Note, that in the case β = β 0 = β 1 = β 2 we have ∂β(t) ∂β = 1. By adding the time delay, we obtain The adjoint system reads as supplemented by the terminal condition (z S , z E , z I , z R , z D )(T Fit ) = 0. In the case of the time delay we receive Here χ [a,b] (t) denotes the characteristic function of the interval [a, b], i.e. we define χ [a,b] (t) = 1 for t ∈ [a, b] and = 0 otherwise. To solve the optimization problem (6a)-(6b) numerically, we apply the Forward-Backward Sweep method [13] combined with a Quasi-Newton method (BFGS) [15].
In each iteration step the ODEs and DDEs of the state variables and adjoint equations are solved with Runge-Kutta methods before the corresponding gradient and direction of descent can be determined. The algorithm stops as soon as the termination condition J(u k+1 ) -J(u k ) < TOL is fulfilled.
As initial values we use β = β 0 = β 1 = β 2 = 0.3 for the transmission rate. This is justified by the fact that an average Basic Reproduction Number of about R 0 = 3 is assumed and in our basic model we have Epidemiologically, R 0 indicates the number of new infections an infected individual causes during the infectious period in an otherwise susceptible population. For the sake of simplicity, we assume the same starting value for I 0 and E 0 . This corresponds to the value at the first data point of our measurement, i.e. 130 registered infected persons on 1st March. As already mentioned, we assume that for the recovered and deaths at this time R 0 = D 0 = 0 holds. The possible problems with the optimization of δ and μ were already mentioned in the previous section. To increase the probability of generating a global minimum, we use n = 1000 normally distributed start values for both parameters fulfilling δ ∼ N (0.25, 0.25 2 ) and μ ∼ N (0.03, 0.03 2 ) with δ, μ > 0. The algorithm selects the best result of these n data fits. The reason for this is the assumption that the proportion of detected cases is between 1-50% and the lethality below 6%. For the case fatality rate μ we have where R reported stands for the reported recovered at time T Fit . The approach for this estimation can be found in [10]. The smaller the detection rate δ, the lower the upper bound for the case fatality gets. In case of the time delay we choose as initial history for s ∈ [t 0τ , t 0 ] This is justified by the fact that the number of registered cases has increased tenfold during this period and we assume an exponential growth in this time span.

Simulation results
To estimate the unknown parameters u, we match the data reported on a daily basis by Johns Hopkins [5] to our simulation results for a time period starting on 1 March. The first results in Fig. 5 show a parameter estimation using the basic model (1a)-(1e) and the time period before the onset of any containment measures, i.e. before the closing of schools on 16 March. We fitted the parameters β, δ and μ along with the initial values E 0 and I 0 over the time period 1 March to 16 March. The initial values E 0 and I 0 are also subject to fitting, since the official data does not provide information about the active infections at a given day. The weight ω 2 = 1 to keep the cost functional dominated by the two least square errors. The other weight is chosen as ω 1 = 500 to compensate the significantly smaller value of the least square error in the fatal cases.  Table 1 Optimal parameter values for the three models (1a)-(1e), (3a)-(3e) and (4a)-(4e) obtained from the minimization problem (6a)-(6b) For the given time period of the fit, the model prediction and the observed data are in good accordance. The estimated parameter values are given in Table 1. The detection rate was estimated as δ = 0.37 implying that the true number of infections exceeds the registered cases by a factor 3. The transmission rate β = 0.57 accounts for a doubling time of 2.6 days at the initial, uncontrolled phase of the epidemic in Germany.
In Fig. 6 we show the results obtained with the time-dependent model (3a)-(3e). In this case, the fitting period equals to the entire simulation period starting from 1 March to 7 April. The weights ω 1 , ω 2 are identical to the previous simulation. The obtained transmission rate and according doubling times change from β 0 = 0.5232 and T 2 (β 0 ) = 2.8 days at the initial uncontrolled phase to β 2 = 0.18 and t 2 (β 2 ) = 11.4 days after the contact ban has been introduced. The effect of the contact ban effectively reduces the transmission rate by a factor of about 3 and significantly slows down the speed of the epidemics by increasing the doubling time by a factor 4.
In Fig. 7 we show the result obtained with the delay model (4a)-(4e). For the delay model, we assume a delay of 14 days between entering the class of infected and death. Again, we show the simulation results compared to the reported cases for the infections and deaths. Quite good agreement is found between the model and the simulation for both, infections and deaths. Compared to the time-dependent model, shown in Fig. 6, the delay model agrees better in particular for the fatal cases. In Table 1 we have listed the estimated parameter values in for the three models. We have also included the normalized L 2 -difference between the simulation outcome and the given data, i.e. the first two summands from the cost fuctional (5). A t-test revealed that the deviations of the simulation to the reported data is not normal distributed at a significance level of 5%.
In the simulations the detection rate is found to be 20-40%, indicating that the true number of SARS-CoV-2 infections might be 3-5 times higher that the officially recorded  Fit of the delay model (4a)-(4e) to the data for the period 1 March to 7 April data suggest. The lethality rate is found to be rather small, taking into account the large number of true cases.
Comparing the obtained values for the lethality, the value for the delay-model seems to be most realistic, since in this model we compare the fatal cases today to the infections that occurred two weeks ago. The two other models related the fatal cases of today to the infected cases today, hence to a significantly larger number. Therefore in these to models, the lethality rate seems to be smaller.

Conclusions and outlook
We present three SIR-based models for describing the outbreak of the SARS-CoV-2 outbreak in Germany. Besides a standard SEIR-model, we consider an extension taking into account the effect of social distancing by a time-dependent reduction of the transmission rate. The third model introduces a delay-term to accurately describe the deaths depending on infected cases that occurred several days in the past. Comparing the simulation results to the data published by Johns Hopkins University allows an estimation of the unknown model parameters. Best results are obtained using the delay equation model. In this setting, we find a detection rate of about 20% and a lethality of about 4%. The social distancing measures were leading to an effective reduction of the transmission rate by a factor 4. That is, after the introduction of the measures roughly just 25% of the social contact compared to the initial period were leading to infections.