Memory-based meso-scale modeling of Covid-19

The COVID-19 pandemic has led to an unprecedented world-wide effort to gather data, model, and understand the viral spread. Entire societies and economies are desperate to recover and get back to normality. However, to this end accurate models are of essence that capture both the viral spread and the courses of disease in space and time at reasonable resolution. Here, we combine a spatially resolved county-level infection model for Germany with a memory-based integro-differential approach capable of directly including medical data on the course of disease, which is not possible when using traditional SIR-type models. We calibrate our model with data on cumulative detected infections and deaths from the Robert-Koch Institute and demonstrate how the model can be used to obtain county- or even city-level estimates on the number of new infections, hospitality rates and demands on intensive care units. We believe that the present work may help guide decision makers to locally fine-tune their expedient response to potential new outbreaks in the near future.


Introduction
The COVID-19 pandemic continues to hold our way of life on this planet in a tight grip. Over the whole world, we have now reached more than 10 million infections [1]. While new infections still rise at alarming pace in the United States, Brazil, or India, most other Asian and European countries that were hit much earlier by the pandemic seem to have succeeded in reducing the number of new daily cases. This success can largely be attributed to fast and locally tailored political measures that introduced severe travel restrictions [2][3][4][5][6] and curtailed public life. Initially, the measures were met by a largely understanding general public. However, the partial necessity of police enforcement and increasing protest against contact restrictions, locally even encouraged by politicians [7], demonstrate rising anger, fear, or even mental health problems due to the current situation [8,9]. Thus, it is critical to carefully reopen the economy and reestablish public life, while avoiding a relapse and a potential collapse of the health-care system, which may entail much stricter measures again. To reach this goal, however, decisions must be made quickly and often locally at county level, based on reliable data and trustworthy predictions. Clearly, accurate models are of essence to capture the disease dynamics at exactly this spatial meso-scale, to predict the number of new infections per day or the number of patients that may require intensive care.
Here, we focus on the situation in Germany, where county-resolved daily and cumulative infection cases are reliably reported by the Robert-Koch Institute (RKI, [10]). We combine two previous modeling advances [11,12] into a locally resolved, history-type model that captures the spatiotemporal evolution of the pandemic in Germany. We use a generalization of typically known SIR-type compartment models that allows for a much better representation of the courses of disease [12]. While becoming infected is well represented by a simple ordinary differential equation (ODE), the remaining course of disease is captured rather restrictively by these ODE-based, SIR-type models [13].
Based on the integro-differential model introduced by Kermack and McKendrick already in 1927 [14] and recently reintroduced by [12,13], we model the spatial spread of Covid-19 in the following way: for all s < 0, t > 0 and x ∈ Ω with Ω ⊂ R 2 open denoting the considered spatial region. The normalized initial history datum is given by ) denotes the normalized number of susceptibles. The weight γ I ∈ L 1 ((0, ∞); R ≥0 )) with γ I L 1 ((0,∞)) = 1 describes the evolution of infectiousness, where γ I (τ ) defines the infectiousness of an individual at τ days after the infection event.
The considered balance law is of nonlocal-history type. Nonlocality as well as history in balance laws are receiving increasing attention to model real world phenomena. They provide a more detailed way to model evolution and can be seen as the mesoscopic link between purely macroscopic and fully microscopic models. In the considered application, the microscopic equivalent-agent models [15][16][17]-can be interpreted as a measure valued solution to the proposed model.
The classically used compartment models (SIR, SEIR) [18] have been widely used to model the viral spread. Their recently revealed relationship to Hamiltonian mechanics is quite insightful [19], demonstrating that they constitute a mere simplification of the here considered integrodifferential equations. In terms of the spatial resolutionwhich of course can also be modeled in the compartment models [20,21]-the classical SIR model can be seen as the singular limit of the interaction term, i.e. β(·, * , ) → b(·)δ( * − ) for a given b ∈ L ∞ ((0, ∞); R ≥0 ). The models are generalized with respect to the evolution of infectiousness of infected individuals. The considered model can represent-based on medical data-any course of infectious-ness, in contrast to, e.g., an assumed exponential decay in the widely used SIR model. As introduced in [11], we discretize our spatial domain Ω, Germany, at county-level (or even city-level), where current containment rules are steadily evaluated and adapted in case local infection numbers rise up again. Our county-interaction network is adapted from the GLobal Epidemic And Mobility (GLEAM) model [22,23], focusing on mid-and shortrange interactions motivated by the severely restricted air travel [24,25]. Taken together, this spatially resolved integrodifferential model allows us to accurately analyze and predict disease dynamics at its various stages and the effect of local measures.

Memory-based spatial infection model
To model the spread of the disease in a discretized spatial setting, we consider a finite partition of the domain Ω, i.e.
where N denotes the number of counties or cities in Germany, depending on the spatial resolution. We obtain the following memory type vector-valued initial value probleṁ is the vector-valued normalized initial history datum. We introduce • to denote the transformation of a vector • into a quadratic diagonal matrix, where the entries along the diagonal equal those of the vector. The time-dependent, vector-valued function S ∈ W 1,∞ ((0, ∞); [0, 1] N ) denotes the normalized number of susceptibles. The matrix-valued function B ∈ L ∞ ((0, ∞); R N ×N ≥0 ) denotes the infection rates and interaction between the considered regions Ω k with k ∈ {1, . . . , N }. The existence and uniqueness of a solution of the proposed integro-differential equations-continuous as well as discretized in space-is proven e.g. in [12] for all γ I for which there exists > 0 s.t. γ I | (0, ) ≡ 0. This is a rather natural condition, since the incubation time-the period during which the infected are not yet infectious-is positive.
Based on the history of S, other quantities and subgroups can be determined directly from including medical data on the various courses and infectiousness levels of the disease via corresponding integration weights: We distinguish between the states infectious γ I , symptomatic γ S , tested and quarantined γ Q , hospitalized γ H , in intensive care γ ICU , recovered γ R and deceased γ D . Following the contribution of [ disease progression in our model: (a) light symptoms, recovering without hospitalization, 95% share; (b) hospitalization, recovering without intensive care, 4% share; (c) patients in intensive care and recovering, 0.4% share; (d) patients in intensive care that eventually die, 0.6% share. Figure 1 depicts the four different courses of the disease as represented in our model, including their corresponding state transitions and infectiousness levels. Note that only the weighted sum γ I = 4 i=1 share i · γ I ,i is necessary for the solution of Eq. (1), but the individual contributions allow for detailed descriptions of disease progression from medical data and corresponding post processing. We normalize the integral over γ I such that it represents the probability of infection.
We further assume that patients in courses (b) to (d) will be tested positive and are thereby considered in reported infection numbers. The ratio of total versus detected infections is defined as dark ratio ω, thereby representing the factor of unknown cases. Since the dark ratio is not necessarily constant in space, this is taken into account by locally scaling the function γ Q that represents the detected and quarantined state of course (a) to the appropriate value. Since the dark ratio seems to closely correlate with testing capacities [11], we introduce federal-state-wise dark ratios ω j , assembled in the vector ω, that vary only over states with j ∈ {1, . . . , 16}, Due to locally differing behavior patterns and, in particular, political measures to reduce social contacts, the infection rates vary in space and time. Based on our previous findings in [11], we introduce federal-state-wise initial infection rates β j0 with j ∈ {1, . . . , 16}, and two reduction factors β red 1 and β red 2 representing the major restrictions of 1) cancelling large events and 2) contact restrictions. Those model parameters are calibrated using data reported by the RKI, as described in Sect. 2.3. Since the shut-down measures were introduced at slightly different times T j1 and T j2 in the different federal states of Germany, we model the time-dependent reduction of infection rates via piece-wise constant functionsβ j1 (t) andβ j2 (t) to obtain the overall infection rates in each state (2) To model cross-county interactions, we adapt the GLEAM short-and mid-range mobility network [22,23] as introduced in [11] and capture cross-county infections by where β k are time-dependent infection rates, c k are the crosscounty infection weights, N k is the number of inhabitants in the largest city of county k, N max = 3e6 corresponds to the number of inhabitants in Germany's largest city Berlin, and r kl is the distance between counties k and l. Importantly, β k and c k are identical for all counties within one federal state, such that Eq. 3 introduces 16 additional model parameters c j with j ∈ {1, . . . , 16} which need to be calibrated based on reported data in the literature (see Sect. 2.3). The parameters for the GLEAM model are taken from [22] and given in Table 1. Figure 2b displays the mobility network across Germany [11]. Note that both the county-internal as well as the cross-county infections contribute to the basic reproduction number in our spatial model. For the distinction of individual counties, we use the municipal directory (Gemeindeverzeichnis) from the German Federal Statistical Office (Statistisches Bundesamt) [27], which delivers area, shape and population data. We consider city-wise population data, which is accumulated over the entire county or corresponding spatial domain for the respective model. The center of population of each county serves as its spatial coordinates.
The detailed description of the courses of disease allows for elaborate post-processing of the solution to evaluate any described quantity. Generally, evaluation differs for cumulative and current quantities. The number of cumulative discovered infections Q or the number of deceased D are evaluated by double integrals such as Current quantities like infectious I(t), those with symptoms, hospitalized or ICU patients follow from expressions such as

Initial conditions
From data the number of positively tested people on day zero is known, but the integro-differential equation model requires initial values for the infected at each spatial node as well as an initial history as a starting point for the integration. Initially, we assume exponential growth in all counties described by the ansatz with ν = 0.345 from fitting an exponential function to RKI data on cumulative COVID-19 cases in all of Germany from March 2 to March 6. The initial estimated number of infected ε at time t = 0 in each county can be calculated by combining Eq. (4), the number of initially reported cases Q 0 , Eq. (6) and the time derivativeṠ. From the result and Eq. (6), the initial history can be estimated.
The high-resolution network model brings with it the challenge for spatially consistent initial conditions. However, most counties did not yet have any known cases on the very first day, limiting the possibility of simply scaling overall initial infections per county by the dark ratio. Thus, for the county-based model we selected the distribution of initial infections according to data for March 16 [10], linearly scaling down the overall number of infections to the number reported on our starting date March 2.

Fitting procedure
To approach our spatially resolved county-model, we followed a cascade optimization strategy. Data analysis and preliminary simulations had shown that we require federalstate-dependent dark ratios ω j and infection rates β j0 , j ∈ {1, . . . , 16} [11]. Figure 2a illustrates the estimates for the state-wise dark ratios ω j , which we obtain by assuming a Germany-wide identical mortality of μ = 6‰ [26] and fitting to the individually reported death tolls, with ω ≈ 6.5. Using state-wise identified dark ratios ω j , we first used a coupled system of 16 nodes connecting each federal state to obtain a Germany-wide average β and reduction factors β red 1 and β red 2 by fitting the cumulative data for Germany via the following objective function The residual R 1 is minimized using a particle swarm optimization (PSO) algorithm described in detail in Sect. 2.4 with weights w I = 1/(t end − t start )/max(Q RKI ), w D = 1/max(D RKI ) and w S = 0.1/max(Q RKI ). We then considered state-wise data to fit β j , j ∈ {1, . . . , 16}, while keeping c = 1, leading to the distribution over Ger-many depicted in Fig. 2c. We fit the cumulative number of confirmed infections reported by the RKI [10] for the time period from March 2 until April 25 with the cumulative number of detected infections Q as defined in Eq. (4), normalized by the maximum cumulative number of reported infections from the RKI. This is the time period during which the various shutdown measures were in place without any noticeable relaxation. On top of that, we include the change-rate of infections on our last day April 25 into the residual vector with the weights The weights are state-wise normalized to balance the contribution of heavily and less affected states. The residual R 2 is again minimized using the PSO algorithm. Finally, we increased the resolution to full county level, amounting to a coupled system of 401 nodes. To re-balance the changed influence of the larger network, we iteratively fitted statewise cross-county weights c j for ∈ {1, . . . , 16} to match the state-wise cumulative infections of the 16 node state-wise model (Q sw ) with the accumulated numbers from the 401 node county-based model (Q cw ) on the last day of the fit. We used a damped gradient-descent like algorithm to update c j at iteration i + 1 following the rule Empirically, we obtained converged cross-county weights within 30 iterations with a limited step size Δc max = 0.25 and a damping exponent ζ = 1.5. The final state-wise distribution of optimized cross-county weights c j is displayed in Fig. 2d.

Particle swarm optimization
Particle swarms are distributed optimization schemes that treat each realization of the d optimization variables as particles with a position x i and a velocity v i in a d-dimensional bounded search space. Particles are initialized with a uniformly random position within the boundaries of the search space and zero initial speed. For the iteration i > 0 the following set of equations describes the behaviour of any particle: The velocity v i+1 is a linear combination of three quantities. The previous velocity v i weighted with the constant factor a results in an inert motion property. The term p i loc − x i represents a force that pulls the particle towards its local attractor p i loc , which is the best position this specific particle has visited so far. Multiplication with the constant weight b loc controls the influence of this quantity. In addition to that the randomized diagonal matrix R i loc with values between 0 and 1 enables optimization in varying directions. The global attractor p i glob represents the best position any particle has visited so far and works analogously to the local attractor with the factor b glob .
We chose established values for a, b loc and b glob as summarized in Table 2, used a total of 300 particles and followed the 'nearest' strategy when particles cross boundaries of the search space during optimization [29]. To prevent overly fast convergence to a visited attractor without broad coverage of the search space, we employed a so-called ring topology neighborhood, such that the global attractor of a particle corresponds only to the best local attractor of its two neighbors below and above. This way, good positions are slowly propagated through the whole swarm, allowing for enhanced exploration of the search space, which well balances run-time efficiency and identification of the true global optimum.

Statistical analysis
To validate the model, we evaluated the temporal correlation between model predictions and RKI data by computing the Pearson correlation coefficient r P , the coefficient of determination R 2 = r 2 P and the corresponding p-value to assess statistical significance via the function [r P , p] =   Figure 3 shows how the optimized spatially resolved memorybased model with 401 network nodes representing each county of Germany well reproduces the cumulative confirmed cases in each of its federal states from March 2 until April 25. For cumulative infection data reported by the RKI [10], we find astonishing and statistically significant (all p < 1e − 12) agreement on the temporal evolution. The only state with an R 2 < 0.98 is Bremen-a city-state with overall very low infection numbers and a population of less than 700.000. Here our quasi-continuum modeling approach and the underlying exponential growth seem to approach their validity limit, and stochastic effects start to prevail. Although only the last data point of reported deaths was considered for parameter identification, the model captures the temporal evolution of Covid-19 related deaths in each state of Germany with remarkable accuracy (all R 2 > 0.91). Here, we observe least agreement in the city-state Hamburg. In general, the model better captures the evolution in higher-populated states, with overall more infections and death tolls. We note that our fitting procedure only operates on state-based information. To further validate our model, we compare county-wise cumulative infection numbers Q as reported by the RKI and our model (Fig. 4 left and Figure 5 shows how the model informs on the temporal evolution of cumulative confirmed cases, with more detailed resolution on the subgroups of symptomatic, infectious, and hospitalized, patients in the ICU, as well as the dead. A first kink in the infectious group is clearly visible at the beginning of March due to the cancellation of major events, which then drops significantly when contact restrictions become effective shortly afterwards. Figure 6 shows the model predicted spatial distribution at county resolution of infectious, symptomatic, hospitalized, and patients in intensive care, following from the individual disease courses in Fig. 1. We consider a period from early March, where the exponential growth of the disease started in Germany, until early June under the assumption that the contact reduction factors stay in place. In early March, most of the infected were at an early stage of the disease, i.e., most of them were infectious but did not have disease specific symptoms yet (on average, the first symptoms appear on the fifth day after the infection event [26]). This explains the delay in symptomatic infections clearly visible in Fig. 5.

Model predictions
In our model, we assume that most of the symptomatic voluntarily quarantine themselves and no longer infect others, implying that the infectiousness decreases when people move to the symptomatic group (cf. Fig. 1). The infectious state ends at the latest when the symptomatic have been tested positive for the virus and are quarantined. The symptomatic state of Covid-19 lasts approximately nine days on average [26], explaining why the symptomatic group is about double in size compared to the infectious group in Fig. 5. Figure 6 also shows the delay in Covid-19 cases that need inpatient treatment or even intensive care. As reported in [26], infected are typically hospitalized for nine days after the infection event at a probability of 4.5%. As this is encoded in the courses of disease, the snapshot on March 2 reports hardly any hospitalized patients. According to [26], patients  Finally, we show how our model can be adapted to locally increase resolution to individual city-or community level. Fig. 7 Spatial distribution of the infectious (I) on April 2 at county level for all of Germany (top) and with locally increased resolution to community-level (bottom). The non-densified part of the domain is greyed out for the sake of better visual contrast. Zoomed regions show county-and community resolution for counties Erlangen, Fürth, Nürnberg and their rural surroundings. Note that the proposed macroscopic model reaches its validity limit for very low daily new infections within one subregion Figure 7 shows the Germany-wide county-level simulation (top), with a zoom into the metropolitan area of Nürnberg, Erlangen and Fürth and its surrounding counties. Increasing the resolution within this domain to community level (bottom) but maintaining county-level for the rest of Germany leads to a network of 464 nodes. The zoom-in clearly shows that county infections are dominated by their largest cities, following the three purple areas that represent Nürnberg, Fürth and Erlangen from bottom to top, underpinning our formulation of the cross-county infection terms (cf. Eq. (3)). Surrounding communities suffer much less infections due to their much smaller populations. Gray areas correspond to rural public space not assigned to a specific community [27].

Conclusions and outlook
We have presented a memory-based network model to predict the spatio-temporal outbreak dynamics of the Covid-19 pandemic in Germany. The model considers the effects of political measures, the cancellation of major events and contact restrictions, and the different possible courses of the disease, which is not possible when using traditional SIRtype models. It well represents the evolution of confirmed cases and deaths reported by the RKI from March 2 until April 25. We have then used the model to predict the further developments until June and have provided estimates for the county-wise required capacity of the local health care system, i.e. the number of patients that require hospitalization and even intensive care. Finally, we have demonstrated that the model can be refined to predict the interaction and local outbreak dynamics at community level.
By now, medical data on observed disease progression at most stages during a COVID-19 infection is abundantly available and continues to improve. Our versatile integrodifferential approach directly integrates these data into the model and can easily be extended, corroborating its superiority over standard SIR-type models. In general, the model can thus handle an arbitrary number of courses of the disease. Similarly, it may expand to consider region-dependent demographics or varying capacities and quality of treatment of the health-care system.
While the model can serve as a valuable tool to assess the effects of new super spreader events-which may occur any time-on the distribution of cases in Germany, it reaches its validity limit when the number of infections becomes small. To additionally capture this even smaller scale, a coupling to individual agent-based models [15][16][17] may be beneficial.