Modelling the Spatial Spread of COVID-19 in a German District using a Diffusion Model

In this study, we present an integro-differential model to simulate the local spread of infections. The model incorporates a standard susceptible-infected-recovered (\textit{SIR}-) model enhanced by an integral kernel, allowing for non-homogeneous mixing between susceptibles and infectives. We define requirements for the kernel function and derive analytical results for both the \textit{SIR}- and a reduced susceptible-infected-susceptible (\textit{SIS}-) model, especially the uniqueness of solutions. In order to optimize the balance between disease containment and the social and political costs associated with lockdown measures, we set up requirements for the implementation of control functions, and show examples for continuous and time-dependent, continuous and space- and time-dependent, and piecewise constant space- and time-dependent controls. Latter represent reality more closely as the control cannot be updated for every time and location. We found the optimal control values for all of those setups, which are by nature best for a continuous and space-and time dependent control, yet found reasonable results for the discrete setting as well. To validate the numerical results of the integro-differential model, we compare them to an established agent-based model that incorporates social and other microscopical factors more accurately and thus acts as a benchmark for the validity of the integro-differential approach. A close match between the results of both models validates the integro-differential model as an efficient macroscopic proxy. Since computing an optimal control strategy for agent-based models is computationally very expensive, yet comparatively cheap for the integro-differential model, using the proxy model might have interesting implications for future research.


Introduction
At the beginning of January 2020, the COVID-19 virus began to spread throughout mainland China, with the consequences that we all have experienced in the last three years.Initially, the number of cases was limited to single clusters in a limited number of locations, but later on expanded throughout the country.In previous studies, we have investigated the macroscopic impact of the epidemic using an SIR-model for all cases in Germany.For this, we have used classical differential models such as the SEIRD-(susceptible-exposed-infected-recovered-dead) models to describe the spread of infections during the first wave (cf.Heidrich et al. 1 ), as well as the impact of travelers on disease dynamics in summer 2020 (cf.Schäfer et al. 2 ), both with a strong emphasis on parameter estimation.
In this study, we aim to model the local spread of infections using PDE (partial differential equation) models to gain a better understanding of the diffusion on smaller scales, similar to the work of Viguerie et al. 3 and Wang and Yamamoto 4 .Kuehn and Mölter 5 also investigated the impact of transportation on epidemics using two coupled models, a static epidemic network and a dynamic transportation network, with non-local, fractional transport dynamics.Logeshwari et al. 6 also provide a fractional PDE model of the spatial spread of COVID-19.A general challenge with diffusion-based PDE models is that diffusion of all compartments leads to unwanted diffusion in the total population, which we aim to avoid.This study uses data down to the municipality level for parameter estimation of the model.
We focus our numerical problem on a single district, the district of Birkenfeld in southwestern Germany within the state of Rhineland-Palatinate. Approximately 81,000 people live there in an area of about 780 km 2 .The largest city within the district is Idar-Oberstein, with about 28,000 inhabitants in an area of about 92 km 2 .The remaining people live in the municipalities of Birkenfeld, Baumholder, and Herrstein-Rhaunen.Within a 1.5 hour drive via federal highways and freeways, the following metropolitan areas can be reached: Mainz, Trier, Koblenz, Kaiserslautern, Saarbrücken, and Frankfurt.In addition, the Frankfurt-Hahn airport is located in the neighboring Rhein-Hunsrück district to the northeast.The region is very rural, with daily commuter movements common in the direction of the aforementioned metropolitan areas.The region is also visited by tourists due to the gemstone industry and trade in Idar-Oberstein, the numerous hiking routes, and the nearby Hunsrück-Hochwald National Park.A map of this district can be found in Fig. 1.
The first COVID-19 case in the district of Birkenfeld was confirmed on March 16, 2020.Until June 30, 2020, there were only 90 registered cases in the entire district.However, the number of cases increased during the second wave in autumn/winter 2020/21, with a cumulative 2,513 cases confirmed until March 31, 2021 7 .As of October 2022, over 32,000 cases have been counted in the district.For several reasons, we restrict our research to the data from the second wave: there was only a limited number of previous infections, but comparatively high infection numbers further on; also, no vaccines were available until the beginning of 2021, as well as a very small amount of persons having been exposed to the disease multiple times.Finally, the political lockdown restrictions, particularly in November and December 2020, led to a slower mixing of cases between different districts.This inter-district mobility is not taken into account by our models.We consider daily infection data from the district and all of its municipalities in the time frame from October 1, 2020 to February 25, 2021.the cumulative number of infections in the Birkenfeld district is depicted in Fig. 2. In this figure, the first step visible in the data between end of December and beginning of January is mainly caused by a delay of registration of case numbers due to Christmas holidays.The reasons for second one in mid January are unknown and could be related to registration delays within the district.An important fact is that there are only one detected initial infection case on October 1, situated in the city of Idar-Oberstein, so there were no detected cases in the rest of the entire district.
As the pandemic continues to spread globally (as of autumn 2022), the possibility of new variants of the virus or future pandemics highlights the importance of modeling the spread of infections, particularly in their early stages when limited medical data is available.The accuracy of these models depends greatly on the used parameters and the corresponding data.In this study, we present SEIRD-models, which are commonly used in epidemiological simulations, and estimate their parameters using data from the Robert-Koch-Institute (RKI) 7 and private communications with the Birkenfeld district government 8 .We perform the estimations using both adjoint and Metropolis methods and base them on a least square fit between the model output and the reported daily data.

PDE models
To model the spatial COVID spread in the presented areas, we use an epidemiological reaction-diffusion model.For this purpose, we consider a corresponding spatial area (x, y) ∈ Ω in a time period t ∈ T := [0,t end ].We are looking for a function u : V → R m with V = Ω × T , which is twice continuously differentiable on Ω and once continuously differentiable on T , briefly u ∈ C (V, R m ) 2,1 .The following PDE system has to be fulfilled: Here, ∂ t u stands for the component-wise derivative of u in the direction of time, i.e., ∂ t u = (∂ t u 1 , ..., ∂ t u m ) T and ∆ x,y = (∂ xx u 1 + ∂ yy u 1 , ..., ∂ xx u m + ∂ yy u m ) T for the Laplace operator in Ω.The parameter κ describes the diffusivity of the system and the function f (u) contains the epidemiological component(s).As an initial condition, at time t = 0 a function u 0 : Ω → R m is used with u(x, y,t = 0) = u 0 (x, y).In addition, Neumann boundary conditions are used, where ∂ ν u = (∂ ν u 1 , ..., ∂ ν u m ) T stands for the derivative in the direction of the outward pointing unit normal ν and ∂ Ω stands for the boundary of Ω.In terms of context, the latter means that no individual can leave or enter the territory Ω.This seems strange at first, since the district of Birkenfeld in practice can be left or entered by land.On the other hand, we are looking at data sets from a period when profound containment measures had already been taken in the region and social measures, including a significant reduction inter-district mobility, were already implemented.
For epidemiological modelling, we make use of a variant of the SIR model introduced by Kermack and McKendrick 10 , the SEIR-model, and consider compartments as functions S, E, I, R, N ∈ C (V, R) 2,1 , which have the following meanings: • Susceptibles S: Depending on the transmission route, these individuals can become infected with the infectious disease when contact occurs.
• Exposed E: The corresponding indiviuals have already ingested the pathogen, but are not yet infectious because they are still in the latency period.
• Infected I: These individuals are infected with the disease and infectious.Contact with a susceptible individual can therefore lead to transmission of the disease.
• Recovered R: After surviving an infection, individuals are considered recovered.These individuals can no longer transmit the disease or get infected.
• Population N: The total number of individuals.

3/16
For instance, I(x, y,t) indicates the number of infected individuals in the spatial coordinate (x, y) ∈ Ω at time t ∈ T .Based on these presented groups, different epidemiological models can now be derived.We present here the PDE systems of three common models in Tab. 1.The derivation and precise functioning of spatial epidemiological models will not be explained in detail here; for this purpose, reference is made to e.g.Martcheva 11 .At the core of every epidemiological model is the so-called incidence term β (t) N SI, which indicates how many individuals are newly infected with the disease in coordinate (x, y) at time t.The incidence term depends on a time dependent transmission rate β : [0,t end ] → R + .In simple models, this can also be assumed to be a constant parameter, but we assume that the transmission rate may fluctuate over the observed periods due to the stepwise restrictions on the population.The value of the transmission rate β is generally unknown and must be adjusted using the data sets.Another parameter in the models is the recovery rate γ.This is the reciprocal of the time required on average for an individual to recover from the disease.Thus, if we assume that t is in days and an individual takes 10 days to recover, it holds γ = 1 10 .In addition, the SEIR-model contains the parameter ϑ , which is the reciprocal of the latency period, i.e. the time between the uptake of the pathogen into the body and the onset of infectiousness.For example, assuming three days, it holds ϑ = 1 3 .It should be noted here that the latency period does not have to be congruent with the incubation period, as the latter indicates the period of time until the onset of the first symptoms.With regard to COVID-19 in particular, it has been shown that infectivity sets in even before the onset of symptoms (cf.He et al. 12 ).Due to simplicity, we also chose We use the SEIR-model as an example to show how the models are prepared for later data fitting.In the first step, we substitute R = N − S − E − I and thus reduce the system to an SEI-model.It should be noted, however, that also for N a PDE has to be solved for which holds For this reason, we normalize the reduced SEI-model by dividing all rows by N, assuming that the population density mathematically fulfills N(x, y,t) > 0 on V .Defining u 1 := S N , u 2 := E N , u 3 := I N and u := (u 1 , u 2 , u 3 ), we obtain a system as in (1) with Tab. 2 summarizes the results for the presented models.
Table 2. Summary of f (u) for the reduced and normalized models.
To meaningfully include the biological context, it must hold γ, ϑ > 0 and κ ≥ 0 and, as initial condition, u 0 ≥ 0 in Ω.In addition, we assume that there are infected individuals in the area Ω, i.e.Ω I 0 dω > 0. For the reduced and normalized SEImodel, it must then hold Ω u 3,0 dω > 0, using the notation u j (x, y, 0) := u j,0 .The notation N 0 := Ω N(x, y, 0) dω represents the total number of individuals at time t = 0 in Ω.It must be valid that N 0 > 0.Moreover, we define the total population in the area Ω at time t as N : [0,t end ] → (0, +∞) with N (t) := Ω N(x, y,t) dω.Due to the Neumann boundary conditions, we receive using Gauss's Thus, the total population in the domain Ω is constant with respect to time.Analytically, there exists a unique solution for each of the PDE systems (1) with the presented f (u) in Table 2 in conjunction with the mentioned preconditions 13 .Due to the formulation using diffusion for the total population, an equilibrium will only set in when the population density in the entire district is equal.Then, the temporal equilibrium will be analogous to the equilibrium of the system without diffusion.
As already mentioned, certain parameters of the model are known, such as γ and ϑ .The transmission rates β j and the diffusivity κ are usually unknown.For the transmission rate, we assume that the time-dependancy is piecewise constant.Due to 'light' lockdown restrictions from November 2, 2020, and 'stricter' restrictions from December 17, 2020 to the end of the observed time interval, we assume three different time intervals as follows: In addition, due to noisy data sets, the initial conditions u 0 must also be adjusted.For that, we present two approaches to solve this in the following sections.

Crank-Nicholson method for the SEIR-model
For this purpose, we discretize the studied region Ω in x and y directions in equal equidistant step sizes h x and h y , respectively.Also, the time interval T = [0,t end ] is divided into equidistant steps of length τ.In the following we use the notation u n i, j = u(x j , y i ,t n ).The Laplace operator is expressed by finite differences and the Crank-Nicolson scheme reads as The goal is to transform this approach so that u n+1 i, j can be solved with a linear system of equations.The non-linearity of f (u) is solved by evolving f using Taylor expansion for small values τ around the current iteration value u n i, j with f (u n+1 i, j ) = f (u i, j (t n + τ)): If we set r x := κτ/h 2 x and r y := κτ/h 2 y , this leads in eqn.(7) using eqn.(8) to ) is negligible for small values of τ, which leads to the system The system (10) leads to a linear equation system The vectors q n and f n are defined analogously, where l x and l y indicate the number of discretization points in x and y direction with respect to Ω.The square and non singular matrices A and B are defined to contain the Neumann boundary conditions, which are implemented by, e.g.u n+1 k+1, j = u n+1 k, j if u n+1 k, j lies on the boundary of the domain ∂ Ω.The previous refers to the solution of the PDE system of the state variable u.The adjoint system must be solved backwards in time, which leads to the approach where the p(u n i, j ) contains componentwise the corresponding discretized terms of ∂ u j g + ∑ m k=1 z k ∂ u j f k .Proceeding analogously as before yields the system Thus, when solving the linear system of equations, the matrices A and B can also be used, due to the same Neumann boundary conditions.

Finite element method for the SEIR-model
An alternative to the Crank-Nicholson method which is used in the adjoint method, we also present a version of the finite element method which produces similar results.By plugging eqn.(3) in eqn.(1), the PDE system reads as follows: As the diffusion and ODE parts in eqs.( 14) are handled by two different schemes, we can make use of an operator splitting (cf.MacNamara 14 ) to solve the system.E.g., for one time step ∆t, this procedure is as follows: (1) Solve ∂ t u i = κ ∆ x,y u i for t = ∆t 2 with the corresponding initial and boundary conditions for i = 1, 2, 3. (2) Solve ∂ t u i = f i (u) for t = ∆t with the corresponding initial and boundary conditions for i = 1, 2, 3.
(3) Solve ∂ t u i = κ ∆ x,y u i for t = ∆t 2 with the corresponding initial and boundary conditions for i = 1, 2, 3.The equation in ( 2) is a simple ODE equation which can be solved by any standard solver, e.g. the Euler method or the method of Runge-Kutta.To solve the equation in ( 1) and (3) on the domain Ω with initial and homogeneous Neumann boundary conditions on ∂ Ω, we consider its weak form gained by multiplication with a test function v ∈ H 1 0 (Ω); this means, that v vanishes at the boundary, i.e., v ≡ 0 at ∂ Ω.The weak form reads as follows:Instead of u ∈ C 2,1 , we now aim to find a function This infinite-dimensional problem has to be solved numerically by discretization.We aim to find a solution u i,h in a finitedimensional subspace V h solving We define the subspace V h on the chosen grid and linearly independent basis functions ϕ j piecewise over subregions where for (x, y) ∈ Ω k ; otherwise, those functions vanish, i.e., ϕ j (x, y) ≡ 0 for (x, y) ̸ ∈ Ω k , j = 1, 2, 3, 4. Then the weak form a(u i,h , v h ) = 0 reads as follows: and, due to the linearity of a, Then the stiffness matrices A and B are defined by where n represents the row corresponding to ( j * , k * ) and m the column corresponding to ( j, k), which depends on the chosen order within the matrices.More information about this can e.g.be found in 15 .The linear equation system with a mass matrix can be solved by any scheme; e.g., a 4-step Runge-Kutta scheme (cf.e.g. 16).

Discretization of the domain
To define the domain Ω for the partial differential equation model, we consider the geographical data for the district (D), the association communities (AC) and the municipalities (M).Using the free online data filtering tool Overpass Turbo 17 from OpenStreetMap, we extracted the relevant geographical data and created the relevant matrices assuming the map segment as rectangular, which appears reasonable due to the very small size of the investigated window.The relevant data for the starting values of the locations are then equally distributed across all relevant subdomains of Ω.The size of the window is L x × L y = 39.23 km × 56.05 km.Using a step size of h x = L x 100 and h y = L y 100 , the discretization of the area yields the following 101x101 matrices.

Target function
In the following, we present the analysis for the reduced SEI-model.The derivation for the other models is analogous.Furthermore, in order to avoid confusion, we rename (u S , u E , u I ) := (u 1 , u 2 , u 3 ).The objective of this section is to present two methods for data fitting of the presented models to the dataset of Birkenfeld County in Germany.Therefore, we introduce an objective function J : R d × C (Ω, R m ) 2,1 defined by On the one hand we make use of the L 2 -norm with respect to V and Ω, defined for example by ∥g∥ L 2 X = X g 2 dx 1/2 , on the other hand we use the Euclidean norm . The L 2 V -norm in J involves fitting the model u to the respective data sets, whereby u data I : V → R m is an interpolation through the data points.Interpolation is performed linearly with respect to the time axis for each grid point.We decided to fit the model to the data point by point in the L 2 V norm.The reason for this is the necessary analysis in the adjoint method in subsection 3.5.
The Euclidean norm in the objective function corresponds to a regularization term and contains the parameters of the model u, which have to be fitted to the data.For example, if we assume an unknown constant transmission rate β (t) ≡ β and a diffusivity κ, we have χ = (β , κ, δ ).The L 2 Ω -norms are to match the initial conditions of u to the corresponding initial guess u data j,0 , which is in the following simulations assumed to be 0. The same applies to the initial guess of χ.In addition, the objective function includes weights w j with j = 0, 1, 2, whose choice will be explained later.The goal is now to minimize J while satisfying the model constraints, i.e. min χ,u 0 J , subject to PDE system (1) . (24)

9/16
Note, that we only use informations concerning the daily reported new infected individuals.Thus, the data include the percentage of daily new infected individuals which shall be fitted to the incidence term β (t)u S u I in the reduced and normalized SEI-model.In addition, one must assume that only a fraction of those actually infected are reported.Therefore, we introduce a detection rate δ , which is unknown.

Models, parameter bounds and initial values
Allowing the following respective constraints of the fitted parameters: the starting values at time t 0 for the detected cumulated infected u 3 0 can be taken from the statistics, while we assume no initial recovered person u 4 0 = 0 and, for the initial amount of exposed persons, u 2 0 = u 3 0 /2 with similar reasoning as in 1 .The initial number of infected is then defined as Table 3. Orders of magnitude of the initial values for adapting the model to the available data. param.

Metropolis algorithm
The first presented method makes use of a Metropolis algorithm (cf.Metropolis et al. 18 , Gelman et al. 19 or Gilks et al. 20 ) for estimation of parameters in the PDE system ( 14) according to the procedure described in Schäfer and Götz 21 and Heidrich, Schäfer et al. 1 .Using the parameter set u 0 as of Tab. 3 as starting conditions, we assign random draws u new from a normally distributed (and thus symmetric) proposal function q, i.e. u new ∼ q(u new |u i−1 ), in every iteration i.
Using the previously defined Ĵ(u) as of eqn.( 23) as the target distribution, we calculate the approximative distribution by whereby c is an arbitrary real value.For the acceptance probability, it follows In eqn.(27), we can see that the value of c is redundant, as it cancels out in the division.If the sample is accepted with the probability p, we set ûi = ûnew ; with the probability 1 − p, the sample is declined, meaning û = ûi−1 according to Rusatsi 22 or Schäfer and Götz 21 .For parameter estimation using the Metropolis algorithm, we use algorithm 1.

Parameter estimation via adjoint functions
For parameter estimation via adjoint functions z ∈ C (V, R m ) 2,1 we use them in conjunction with a Lagrange function defined as At a possible minimum (χ * , u 0 * , u * , z * ) of problem (24) must apply The derivatives in directions representing functions are determined with the help of Gâteaux derivatives.The application of the optimal control theory and Pontryagin's maximum (minimum) principle provide the following optimality conditions: (ii) u j,0 = u data j,0 − z j (x,y,0) It should be noted that the present conditions are given as g for the reduced SEI-model.The adjoint equations PDE system has be solved numerically backward in time.Let us assume a time-dependent transmission rate of the form (5).This leads to χ = (β 0 , β 1 , β 2 , κ, δ ) and to a corresponding gradient where As adjoint equations we get in this case For parameter estimation, we use algorithm 2.

11/16
Algorithm 2 Pseudocode for the parameter estimation via adjoint functions.← load initial values and data 2: u, z ← solve PDE for state variable and adjoint function 3: J, ∇J ← compute objective function and gradient regarding χ = (β 0 , ..., β k , κ, δ ) 4: s 1 ← compute search direction regarding χ (Quasi-Newton (BFGS)) 5: s 2 ← ( ũ0 − u 0 ) compute search direction for u 0 with ũ j,0 = u data j,0 − z j (x,y,0) w 2 6: repeat 7: J old ← J Regarding the optimization of χ, a Quasi-Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS) search direction is used.The initial condition u 0 is updated with a convex combination between old value and current 'optimal' value.The step size is mediated by the Armijo stepsize rule.In each optimization step, the PDE system for the state u and adjoint z variable must be solved.This is done by the Crank-Nicholson method presented above.The fact that the state variable must be solved forward and the adjoint variable backward in time also leads to the term "forward-backward sweep method".

Without penalty term (Metropolis)
Using the Metropolis algorithm as of chapter 3.4, we firstly set w 0 = 1 and w 1 = w 2 = 0 in eqn.(23).In Fig. 5, the results for the district of Birkenfeld and in Fig. 6, the results for the four ACs are presented; the red dots represent the respective data, the blue line the outcome of the model.

With penalty term
We now include a penalty term; in eqn.(23) we therefore choose w 0 = 1 and w 1 = w 2 = 1 • 10 −5 , which guarantees a convex problem, cf.Heidrich, Schäfer et al. 1 .In Fig. 7, the results for the district of Birkenfeld and in Fig. 8, the results for the four associated communities are presented; the red dots represent the respective data, the blue line the outcome of the model.We also compare the results to those of a standard SEIR-model with the usual parameter estimation (using the Metropolis algorithm).

Comparison
In the results of the parameter estimation as shown in Tab. 4, it is notable that across all we find β 1 ≈ β 2 , which indicates that the more severe lockdown restrictions from December 17 onwards might be overlaid by the rising number of festivities during Christmas.The adjustment regarding the detection rate shows that, according to the model, the actual number of infected people is approximately three to five times higher than registered, which is in line with previous findings.The parameter values for the transmission parameters β i are relatively consistent across the different methods and parameter settings.The values β 0 > β 1 ≈ β 2 show that the so-called 'light' lockdown in November 2020 had the most significant effect on the case values, while the more 'severe' lockdown before Christmas did not have a significant effect.A reason for this can be Christmas itself during which the amount of contacts has rised by nature, it also has to be noted that the measurement of case numbers during or after Christmas was not consistent, which might lead to delays in the infection data.
The parameter κ is ranging significantly between 0.05 and 0.12, but correlations to the higher initial values of the infected, as well as the detection rate δ , vary across the different methods and parameter settings, depending not only on the weights w 1 and w 2 but also the chosen κ.The results indicate that there is some potential variation in the optimal parameter values and parameter cross-correlation of parameters, e.g. of the initial infection values and the detection rate, as well as to some extend in the diffusivity.This can also be seen in the target function J(u) which, despite the deviations, shows only minor variations among the different methods and parameter settings.
The estimations for some of the districts is generally quite accurate, e.g. for Herrstein-Rhaunen and Baumholder.However, for the city of Idar-Oberstein, the area with most infections, the model underestimates the infection cases, which can be explained by the diffusivity of the model that neglects the urban structures.On the other side, the model overestimate the infection cases in Birkenfeld.In total, the amount of infected is slightly underestimated towards the end, while the Metropolis algorithm with c = 10 −5 provides the most accurate results for the cumulative data.Compared to the SEIR-model without diffusion, we see that both Metropolis and adjoint methods provide better target function values, as the SEIR-model generally overestimates the infections in all districts except the largest (Idar-Oberstein).

Discussion
This study presented a reaction-diffusion model used to simulate the spatial spread of COVID-19, making use data down to the municipality level for parameter estimation.The SEIR-model is based on a set of PDEs that describe the dynamics of susceptible, exposed, infected, and recovered individuals, as well as an incidence term that represents the transmission of the disease.The optimal control is based on a least square fit between the model output and the reported daily data.Two different approaches for the estimation of parameters and approximation of the infection data -the Metropolis algorithm and the adjoint method -were described and implemented, and their results were plotted and compared.Regarding the graphical and numerical results, all routines have provided meaningful results.The models depict the infection values quite accurately in several subdistricts, yet slightly over-or underestimate them in others, which can partially be explained by non-homogeneous behaviour of cities compared to rural areas.On a local level, the quality of the estimations decreases -as expected, as it cannot be assumed that a global model will apply perfectly to the behavior of single villages, especially as some of them had no or less than a handful of detected infected in the observed time interval.The parameter values, including initial infection values, detection rate, and diffusivity, vary across different methods and parameter settings, while the transmission-related parameters remain relatively consistent, and the target function are very similar.Compared to a non-spatial SEIR-model, both Metropolis and adjoint method provide better results with respect to the target function J(u).
It is important to note that the PDE model still provides valuable insights and information and manages to describe the diffusion in the epidemiological situation (very few and locally condensed initial cases in a more or less completely susceptible population).Further refinement of the model or the use of additional data sources may improve its accuracy in the future.Nevertheless, the accuracy of the model on the district and ACs level make it a valuable tool for understanding the spread of 15/16 infections.Not only new variants of COVID-19, but also the possibility of future pandemics underscore the need for accurate modeling on local level.

Figure 1 .Figure 2 .
Figure 1.Map of Birkenfeld County with administrative boundaries of central localities 9 .

Figure 3 .Figure 4 .
Figure 3. Left: Discretizations of the whole district.Right: Discretization of the association communities.

Figure 5 .
Figure 5. Result of the optimization with the Metropolis algorithm for the district of Birkenfeld with w 1 = w 2 = 0.

Figure 6 .
Figure 6.Result of the optimization with the Metropolis algorithm for the lower level administrative units (Verbandsgemeinden) with w 1 = w 2 = 0.

Figure 8 .
Figure 8. Result of the optimization with the various methodes for the lower level administrative units (Verbandsgemeinden) with w 1 = w 2 = 10 −5 .

Table 1 .
Basic examples of epidemiological compartment models with flow chart and PDE system.