# Hybrid PDE-ODE Models for Efficient Simulation of Infection Spread in
Epidemiology

Kristina Maier
Zuse Institute Berlin
Berlin, 14195
<EMAIL_ADDRESS>
& Martin Weiser
Zuse Institute Berlin
Berlin, 14195
<EMAIL_ADDRESS>
& Tim Conrad
Zuse Institute Berlin
Berlin, 14195
<EMAIL_ADDRESS>

###### Abstract

This paper introduces a novel hybrid mathematical modeling approach that
effectively couples Partial Differential Equations (PDEs) with Ordinary
Differential Equations (ODEs), exemplified through the simulation of
epidemiological processes. The hybrid model aims to integrate the spatially
detailed representation of disease dynamics provided by PDEs with the
computational efficiency of ODEs. In the presented epidemiological use-case,
this integration allows for the rapid assessment of public health
interventions and the potential impact of infectious diseases across large
populations. We discuss the theoretical formulation of the hybrid PDE-ODE
model, including the governing equations and boundary conditions. The model’s
capabilities are demonstrated through detailed simulations of disease spread
in synthetic environments and real-world scenarios, specifically focusing on
the regions of Lombardy, Italy, and Berlin, Germany. Results indicate that the
hybrid model achieves a balance between computational speed and accuracy,
making it a valuable tool for policymakers in real-time decision-making and
scenario analysis in epidemiology and potentially in other fields requiring
similar modeling approaches.

The source code is available at this URL: https://zenodo.org/records/11121452.

_Keywords_ Partial differential equation $\cdot$ ordinary differential
equation $\cdot$ SEIR model $\cdot$ epidemic modeling $\cdot$ spatial modeling
$\cdot$ coupling approach $\cdot$ diffusion $\cdot$ agent-based model $\cdot$
parameter optimization $\cdot$ Levenberg-Marquardt algorithm

##  1 Introduction

In the rapidly evolving field of mathematical modeling, the ability to
accurately and efficiently simulate large complex systems is of great
importance. This is particularly true in the field of epidemiology, where
understanding the dynamics of disease spread can have significant implications
for public health and policy decisions. In this field, multiple modeling
approaches exist, such as agent-based models (ABMs) [25], models based on
partial differential equations (PDEs) [13], and ordinary differential equation
(ODE) models [36].

While ABMs are well-suited for simulating detailed actions and interactions at
the individual level, they can be computationally very demanding, especially
when simulating large populations over extended periods or aiming at policy
optimization [26]. This reduces their suitability for real-time analysis or
situations in which fast computational processing is essential.

PDE models, on the other hand, provide a complementary perspective to the
high-resolution perspective provided by ABMs. PDE models are generally less
computationally demanding compared to ABMs and thus allow simulations of
larger geographical areas, in a spatially continuous environment. They provide
valuable insights into the general transmission patterns by considering
factors like population density and mobility, which enable modelers to gain
valuable information about the modeled systems on a regional or national
level, such as possible areas of high infection rates. However, solving large
PDE systems also has some considerable computational complexity – which can
become a significant limitation, particularly when dealing with large
environments and extensive scenarios on a large scale.

ODE models, in contrast, emphasize the temporal dynamics of disease
transmission, offering computational efficiency but lacking the spatial
resolution of PDE models and the detailed granularity of ABMs. This makes them
extremely valuable in situations when rapid analysis is needed. Nevertheless,
this benefit is accompanied by the drawback of decreased accuracy in
representing spatially diverse situations.

While each of these modeling approaches — ABMs, PDEs, and ODEs — plays a
crucial role in understanding epidemiological dynamics, integrating these
models to leverage their distinct advantages remains challenging. In this
paper, we introduce a hybrid PDE-ODE model that seeks to merge these
approaches efficiently, providing a new perspective for enhanced
epidemiological modeling.

Coupled PDE-ODE have been used to model various phenomena, such as tumor
development[12], traffic control[23] and others, alongside specialized solvers
[11]. However, not many applications are known in the context of disease
dynamics. In our context, the innovative aspect of this coupling is the
incorporation of the ODE model as a space-invariant component within the PDE
framework. This integration allows for the exchange of boundary information
between the two models at each time step, guaranteeing a dynamic yet
consistent interplay between local and global disease dynamics. The resulting
hybrid PDE-ODE model offers a faster alternative to a full-PDE model while
maintaining an acceptably low loss of accuracy.

In summary, this paper introduces a new approach for coupling partial
differential equation (PDE) and ordinary differential equation (ODE) models,
in this case tailored for epidemiological applications but with potential
implications beyond this field. By achieving a balance between the detailed
spatial representation of PDEs and the computational efficiency of ODEs, this
approach represents a significant advancement in the field of mathematical
modeling. In the field of epidemiology, this provides a powerful tool for
policymakers to rapidly and accurately assess the impact of various
interventions and scenarios in the management of infectious diseases.

The source code is available at this URL: https://zenodo.org/records/11121452.

##  2 Full-PDE Model

In this section, we provide a comprehensive formulation of the hybrid PDE-ODE
model, designed for an epidemiological use-case, in particular for the study
of infection spreading. We begin by detailing the high-fidelity PDE model
(Sec. 2.1) that captures the spatial dynamics of the disease, including the
definition of initial values (Sec. 2.2) to accurately represent the starting
conditions of the disease spread. We then proceed to discuss the approximation
of regions of lesser interest using ODEs (Sec. 2.3), and finally, present the
methods for seamlessly coupling these models. This hybrid approach not only
retains the spatial resolution critical for accurate simulation but also
leverages the computational simplicity of ODEs, resulting in an efficient yet
powerful tool for analyzing epidemiological data. This structured approach
ensures that the model can simulate the progression of infectious diseases
with both spatial detail and computational efficiency.

###  2.1 Model Equations

We consider a spatial spreading model of SEIR type in a domain
$\Omega\subset\mathbb{R}^{2}$. We employ a reaction-diffusion equation
modeling both the SEIR dynamics and spatial diffusion of infectious. The
choice of the SEIR model is motivated by its ability to capture the essential
stages of disease progression (susceptible, exposed, infectious, and removed),
providing a detailed and realistic framework for understanding the dynamics of
disease transmission within a population.

Let $s,e,i,r:\Omega\times[0,T]\to\mathbb{R}$ denote the fraction of population
density of susceptible, exposed, infectious, and removed (recovered and immune
or deceased) subjects, such that $s+e+i+r=1$. Their evolution is given as

$\displaystyle\partial_{t}s$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
s)-\left(1-\frac{A}{n+n_{0}}\right)s(\beta_{e}e+\beta_{i}i)$ in $\Omega$, (1)
$\displaystyle\partial_{t}e$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
e)+\left(1-\frac{A}{n+n_{0}}\right)s(\beta_{e}e+\beta_{i}i)-\sigma
e-\phi_{e}e$ in $\Omega$, $\displaystyle\partial_{t}i$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla i)+\sigma e-\phi_{i}i$
in $\Omega$, $\displaystyle\partial_{t}r$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla r)+\phi_{i}i+\phi_{e}e$
in $\Omega$, $\displaystyle\nu^{T}D\nabla s$ $\displaystyle=\nu^{T}D\nabla
e=\nu^{T}D\nabla i=\nu^{T}D\nabla r=0$ on $\partial\Omega$

with $0<n\in L^{\infty}(\Omega)$ denoting the given population density and
$\nu$ being the unit outer normal to $\Omega$. The diffusion coefficient $D>0$
describes the spatial mixing of the subpopulations and may, in general, depend
on the spatial position. The constants $\beta_{e}>0$ and $\beta_{i}>0$
describe the infection rates and $\phi_{e}$ and $\phi_{i}$ the recovery rates.
The rate of becoming infectious is denoted as $\sigma$. A description of the
parameters can be found in Table 1. The homogeneous Neumann boundary
conditions ensure that the fractions add up to one for all times.

The Allee effect plays a crucial role in (cell) population biology, explaining
the positive correlation between per capita growth rate and (cell) population
density [20, 21]. It can be expressed as $(1-A/n)$, where $A$ represents the
Allee threshold. If the population density $n$ falls below this threshold $A$,
the net growth rate becomes negative, indicating potential extinction. As the
population density approaches the threshold, growth slows down but remains
positive. Once the population density surpasses the threshold, the Allee
effect diminishes, making it practically irrelevant when $A\ll n$. In [33,
34], the Allee effect was instead used to model the infection rate as a
function of population density, where it represents a lower transmission
probability in less densely populated regions. This can simply be referred to
as a nonlinear infection rate [32] and is similar to incidence functions [6].
The influence of population density on the infection rate is controversial [7,
30, 35]. Nevertheless, we adopt the effect and bound it from below by
$\nicefrac{{1}}{{3}}$. This ensures that the effect is at most three times
lower in sparsely populated areas than in regions with high population
density. To enforce the lower bound, we apply a shift $n_{0}$ in the Allee
term. We choose $n_{0}$ such that the factor is bounded from below

$\displaystyle 1-\frac{A}{n+n_{0}}\geq\frac{1}{3}\qquad\Leftrightarrow\qquad
n+n_{0}\geq\frac{3}{2}A$

for all $n\in\mathbb{R}^{+}$. Hence, we choose
$n_{0}\coloneqq\nicefrac{{3}}{{2}}A$.

parameter | definition
---|---
$s(x,t)$ | fraction of population density of susceptibles
$e(x,t)$ | fraction of population density of exposed
$i(x,t)$ | fraction of population density of infectious
$r(x,t)$ | fraction of population density of removed
$n(x,t)$ | population density
$\beta_{i}$ | symptomatic infection rate
$\beta_{e}$ | asymptomatic infection rate
$\sigma$ | inverse of incubation period
$\phi_{e}$ | asymptomatic recovery rate
$\phi_{i}$ | infected recovery rate
$D$ | diffusion coefficient
Table 1: Parameters of the model (1).

###  2.2 Initial Conditions

Next, we address the challenge of defining initial starting values and
estimating distribution patterns for our model. Due to the typical
unavailability of precise spatial densities for exposed, infected, and removed
populations, we use aggregated data, which is usually available at broader
geographical levels such as districts, subregions, or provinces. The following
outlines our approach for setting initial conditions for fractional densities
within the PDE system. This step is important to ensure that our simulations
start with an accurate and realistic representation of the modeled disease’s
spatial distribution – based on the available data.

Let $\Omega=\bigcup_{l=1}^{L}\Omega^{(l)}$ be a decomposition of $\Omega$ into
provinces $\Omega^{(l)}$ for which aggregated numbers of infectious and of
removed are given.

The initial conditions are computed using the real (smoothed) number
$\hat{\mathcal{I}}^{(l)}$ of infectious and $\hat{\mathcal{R}}^{(l)}$ of
removed individuals for each subdomain $\Omega^{(l)}$, $l=1,\dots,L$. Our
given data for each subdomain $\Omega^{(l)}$ is

$\displaystyle\hat{\mathcal{N}}^{(l)}$ $\displaystyle\text{(total population
number)},$ $\displaystyle\hat{\mathcal{I}}^{(l)}$ $\displaystyle\text{(number
of infectious individuals)},$ $\displaystyle\hat{\mathcal{R}}^{(l)}$
$\displaystyle\text{(number of removed individuals)}.$

Further, we set

$\displaystyle\hat{\mathcal{E}}^{(l)}=\Sigma_{\mathcal{E}}\hat{\mathcal{I}}^{(l)}$
$\displaystyle\text{(number of exposed)},$
$\displaystyle\hat{\mathcal{S}}^{(l)}=\hat{\mathcal{N}}_{l}-(1+\Sigma_{\mathcal{E}})\hat{\mathcal{I}}^{(l)}-\hat{\mathcal{R}}^{(l)}$
$\displaystyle\text{(number of susceptibles)},$

where $\Sigma_{\mathcal{E}}$ scales the number of infectious to estimate the
number of exposed individuals. We use isotropic Gaussian functions
$G^{(l)}(x)=\tfrac{1}{2\pi\sigma^{2}}\exp\Bigl{(}-\tfrac{\|x-x_{0}^{(l)}\|^{2}}{2\sigma^{2}}\Bigr{)}$
located at the centers $x_{0}^{(l)}$ of the provinces
$\Omega^{(l)}\subset\Omega$ for representing density and density fractions as

$n(x)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{n}^{(\tilde{l})}G^{(\tilde{l})}(x),\quad
s(x,0)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{s}^{(\tilde{l})}G^{(\tilde{l})}(x),\quad
i(x,0)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{i}^{(\tilde{l})}G^{(\tilde{l})}(x),\quad\text{and}\quad
e(x,0)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{e}^{(\tilde{l})}G^{(\tilde{l})}(x).$

The coefficients $w^{(\tilde{l})}$ are determined such that the integral of
the densities over each subdomain area $\Omega^{(l)}$ is equal to the given
target data
$\hat{\mathcal{N}}_{l},\hat{\mathcal{S}}^{(l)},\hat{\mathcal{I}}^{(l)},\hat{\mathcal{E}}^{(l)}$
at time $t=0$, i.e.,

$\displaystyle\int_{\Omega^{(l)}}n(x)\;dx=\hat{\mathcal{N}}_{l},$
$\displaystyle\int_{\Omega^{(l)}}s(x,0)n(x)\;dx=\hat{\mathcal{N}}_{l}-(1+\Sigma_{\mathcal{E}})\hat{\mathcal{I}}^{(l)}-\hat{\mathcal{R}}^{(l)},$
$\displaystyle\int_{\Omega^{(l)}}i(x,0)n(x)\;dx=\hat{\mathcal{I}}^{(l)},$
$\displaystyle\int_{\Omega^{(l)}}e(x,0)n(x)\;dx=\Sigma_{\mathcal{E}}\hat{\mathcal{I}}^{(l)},$

leading to linear equation systems of dimension $\tilde{L}$. As $n$ is
required for computing the initial densities $s,i$, and $e$, it is computed
first.

###  2.3 Hybrid PDE-ODE Model

In the following section, we describe how we define our hybrid model based on
the PDE model by segmenting the modeling domain into a PDE part and a – less
complex – ODE part. Our approach is inspired by previous work where similar
concepts have been used for models simulating microbial sensing [19] and
airborne pathogen movement in an SIR framework [9]. Our approach adopts these
concepts to simplify computational demands without compromising the essential
dynamics of disease spread. We utilize integro-differential equations for the
coupling, enabling us to maintain high-resolution modeling where it is crucial
while simplifying other areas to enhance computational efficiency. This will
ultimately allow us to efficiently model large-scale epidemiological phenomena
by reducing complexity in areas where high resolution is less critical. We
detail the mathematical underpinnings of this reduction, including the
seamless integration and boundary conditions that ensure continuity and
consistency across the coupled models.

We first consider a non-overlapping domain decomposition [5, 27]. We assume
the domain $\Omega$ comprises a region $\Omega_{1}\subset\mathbb{R}^{2}$ of
interest to be modeled in spatial detail, and an adjacent region
$\Omega_{2}\subset\mathbb{R}^{2}$ of less interest, where a coarser model is
sufficient. In the region of interest, we employ a reaction-diffusion equation
modeling both the SEIR dynamics and spatial diffusion of infectious. On the
domain $\Omega_{2}$, we assume homogeneous mixing and formulate a
straightforward ODE model. The models are coupled across their joint boundary
$\Gamma:=\partial\Omega_{1}\cap\partial\Omega_{2}$. This model coupling aims
to achieve significantly reduced computational costs compared to implementing
the PDE model on the entire domain
$\Omega\coloneqq\mathop{\mathrm{missing}}{int}\overline{\Omega_{1}\cup\Omega_{2}}$.

Let us first consider the structure of the full-PDE solution on $\Omega$.
Across the interface $\Gamma$ separating $\Omega_{1}$ from $\Omega_{2}$,
continuity of densities and continuity of fluxes hold almost everywhere.
Introducing $s_{1}=s|_{\Omega_{1}}$ and $s_{2}=s|_{\Omega_{2}}$, this reads

$\displaystyle s_{1}$ $\displaystyle=s_{2}$ $\displaystyle\text{a.e. on
$\Gamma$},$ (2) $\displaystyle\nu_{1}^{T}D\nabla s_{1}$
$\displaystyle=\nu_{1}^{T}D\nabla s_{2}$ a.e. on $\Gamma$ (3)

with $\nu_{1}$ being the outer unit normal of $\Omega_{1}$. The same holds for
$e$, $i$, and $r$.

Extending $s_{1}$ and $s_{2}$ by zero to the whole domain $\Omega$, we can
write the relative density of susceptibles in $\Omega$ as $s=s_{1}+s_{2}$.
Multiplying the first PDE of model (1) by a test function $\varphi$ and
integrating over $\Omega_{2}$ results in

$\displaystyle\int_{\Omega_{2}}\dot{s}\varphi\;dx$
$\displaystyle=\int_{\Omega_{2}}\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
s)-\left(1-\frac{A}{n+n_{0}}\right)s(\beta_{e}e+\beta_{i}i)\bigr{)}\varphi\;dx$
$\displaystyle\Leftrightarrow\int_{\Omega_{2}}\dot{s}_{2}\varphi\;dx$
$\displaystyle=\int_{\Omega_{2}}\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
s_{2})-\left(1-\frac{A}{n_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})\bigr{)}\varphi\;dx$
$\displaystyle=\int_{\Omega_{2}}-\nabla\varphi^{T}D\nabla
s_{2}-\left(1-\frac{A}{n_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})\varphi\;dx+\int_{\Gamma}\varphi\nu_{2}^{T}D\nabla
s_{2}\;do$ $\displaystyle=\int_{\Omega_{2}}-\nabla\varphi^{T}D\nabla
s_{2}-\left(1-\frac{A}{n_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})\varphi\;dx-\int_{\Gamma}\varphi\nu_{1}^{T}D\nabla
s_{1}\;do.$

Further, let us decompose $s_{2}=\bar{s}_{2}+\tilde{s}_{2}$ into a spatially
constant part $\bar{s}_{2}=\frac{1}{|\Omega_{2}|}\int_{\Omega_{2}}s_{2}\;dx$
and a fluctuating part $\tilde{s}_{2}$ with mean zero.

Testing with $\varphi=\frac{1}{|\Omega_{2}|}$ yields

$\displaystyle\frac{1}{|\Omega_{2}|}\int_{\Omega_{2}}\dot{s}_{2}\;dx$
$\displaystyle=\frac{1}{|\Omega_{2}|}\int_{\Omega_{2}}-\left(1-\frac{A}{\bar{n}_{2}+\tilde{n}_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})\;dx-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do$ $\displaystyle\Leftrightarrow$ $\displaystyle\dot{\bar{s}}_{2}=$
$\displaystyle\frac{1}{|\Omega_{2}|}\biggl{(}\int_{\Omega_{2}}-\bigl{(}1-A(\bar{n}_{2}+\underbrace{\tilde{n}_{2}}_{\mathclap{\text{small
assuming homogeneous
mixing}}}+n_{0})^{-1}\bigr{)}\bigl{(}\bar{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})+\underbrace{\tilde{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})+\bar{s}_{2}(\beta_{e}\tilde{e}_{2}+\beta_{i}\tilde{i}_{2})}_{=0\text{
if integrated}}$
$\displaystyle+\underbrace{\tilde{s}_{2}(\beta_{e}\tilde{e}_{2}+\beta_{i}\tilde{i}_{2})}_{\mathclap{\text{small
assuming homogeneous mixing}}}\bigr{)}\;dx-\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do\biggr{)}$ $\displaystyle=$
$\displaystyle-\bar{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})\frac{1}{|\Omega_{2}|}\int_{\Omega_{2}}\left(1-\frac{A}{\bar{n}_{2}+n_{0}}\right)\;dx-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do$ $\displaystyle=$
$\displaystyle-\left(1-\frac{A}{\bar{n}_{2}+n_{0}}\right)\bar{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do.$

Thus, we obtain an ODE for $\bar{s}_{2}$ with a right hand side contribution
from the flux of the PDE solution $s_{1}$ across $\Gamma$. The same derivation
applies to $i$ and $e$. Consequently, the mean subpopulation density fractions
in $\Omega_{2}$ evolve according to the ODE model

$\displaystyle\dot{\bar{s}}_{2}$
$\displaystyle=-\left(1-\frac{A}{\bar{n}_{2}+n_{0}}\right)\bar{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do,$ $\displaystyle\dot{\bar{e}}_{2}$
$\displaystyle=\left(1-\frac{A}{\bar{n}_{2}+n_{0}}\right)\bar{s}_{2}(\beta_{e}\bar{e}_{2}+\beta_{i}\bar{i}_{2})-\sigma\bar{e}_{2}-\phi_{e}\bar{e}_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
e_{1}\;do,$ $\displaystyle\dot{\bar{i}}_{2}$
$\displaystyle=\sigma\bar{e}_{2}-\phi_{i}\bar{i}_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
i_{1}\;do,$ $\displaystyle\dot{\bar{r}}_{2}$
$\displaystyle=\phi_{i}\bar{i}_{2}+\phi_{e}\bar{e}_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
r_{1}\;do.$

The statement in [17] about the Allee term not being included in the ODE model
due to the lack of spatial variation in population density overlooks an
important aspect. The Allee term has a maximum value of 1 and is typically
smaller, even negative in their case. Therefore, it is crucial to consider the
Allee term as a mean across the area in the ODE model, ensuring that the
infection numbers in the ODE model do not explode and are more comparable to
the results of the PDE model.

Note that we already assumed $\tilde{s}_{2}$, $\tilde{e}_{2}$, $\tilde{i}_{2}$
and $\tilde{r}_{2}$ to be zero as they get vanishingly small due to
homogeneous mixing. Further, assuming that the total number of susceptibles
$S_{2}(0)$, exposed $E_{2}(0)$, infectious $I_{2}(0)$ and removed $R_{2}(0)$
of the ODE model at time $t=0$ are given, and therefore the (constant) total
population size $N_{2}$, we can set the initial conditions as

$\displaystyle\int_{\Omega_{2}}s_{2}(x,0)n_{2}(x)\;dx=S_{2}(0)$
$\displaystyle\Leftrightarrow$
$\displaystyle\int_{\Omega_{2}}\bar{s}_{2}(0)n_{2}(x)\;dx=S_{2}(0)$
$\displaystyle\Leftrightarrow$
$\displaystyle\bar{s}_{2}(0)=\frac{S_{2}(0)}{N_{2}}.$

Turning to $\Omega_{1}$, the density continuity (2) implies the Dirichlet
boundary conditions $s_{1}=s_{2}$ on $\Gamma$ of the hybrid model. For the
sake of completeness, the hybrid model is presented in the following:

$\displaystyle\partial_{t}s_{1}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
s_{1})-\left(1-\frac{A}{n_{1}+n_{0}}\right)s_{1}(\beta_{e}e_{1}+\beta_{i}i_{1})$
in $\Omega_{1}\times]0,T[$, (4) $\displaystyle\partial_{t}e_{1}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
e_{1})+\left(1-\frac{A}{n_{1}+n_{0}}\right)s_{1}(\beta_{e}e_{1}+\beta_{i}i_{1})-\sigma
e_{1}-\phi_{e}e_{1}$ in $\Omega_{1}\times]0,T[$,
$\displaystyle\partial_{t}i_{1}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla i_{1})+\sigma
e_{1}-\phi_{i}i_{1}$ in $\Omega_{1}\times]0,T[$,
$\displaystyle\partial_{t}r_{1}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
r_{1})+\phi_{e}e_{1}+\phi_{i}i_{1}$ in $\Omega_{1}\times]0,T[$,
$\displaystyle\dot{s}_{2}$
$\displaystyle=-\left(1-\frac{A}{n_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
s_{1}\;do$ in $]0,T[$, $\displaystyle\dot{e}_{2}$
$\displaystyle=\left(1-\frac{A}{n_{2}+n_{0}}\right)s_{2}(\beta_{e}e_{2}+\beta_{i}i_{2})-\sigma
e_{2}-\phi_{e}e_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
e_{1}\;do$ in $]0,T[$, $\displaystyle\dot{i}_{2}$ $\displaystyle=\sigma
e_{2}-\phi_{i}i_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
i_{1}\;do$ in $]0,T[$, $\displaystyle\dot{r}_{2}$
$\displaystyle=\phi_{i}i_{2}+\phi_{e}e_{2}-\frac{1}{|\Omega_{2}|}\int_{\Gamma}\nu_{1}^{T}D\nabla
r_{1}\;do$ in $]0,T[$, $\displaystyle s_{1}$ $\displaystyle=s_{2},\qquad
e_{1}=e_{2},\qquad i_{1}=i_{2},\qquad r_{1}=r_{2}$ a.e. on
$\Gamma\times]0,T[$, $\displaystyle\nu_{1}^{T}D\nabla s_{1}$
$\displaystyle=\nu_{1}^{T}D\nabla e_{1}=\nu_{1}^{T}D\nabla
i_{1}=\nu_{1}^{T}D\nabla r_{1}=0$ on
$(\partial\Omega_{1}\setminus\Gamma)\times]0,T[$, $\displaystyle s_{1}(x,0)$
$\displaystyle=\sum_{\tilde{l}=1}^{\tilde{L}}w_{s}^{(\tilde{l})}G^{(\tilde{l})}(x),\qquad
e_{1}(x,0)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{e}^{(\tilde{l})}G^{(\tilde{l})}(x)$
in $\Omega_{1}$, $\displaystyle i_{1}(x,0)$
$\displaystyle=\sum_{\tilde{l}=1}^{\tilde{L}}w_{i}^{(\tilde{l})}G^{(\tilde{l})}(x),\qquad
r_{1}(x,0)=\sum_{\tilde{l}=1}^{\tilde{L}}w_{r}^{(\tilde{l})}G^{(\tilde{l})}(x)$
in $\Omega_{1}$, $\displaystyle s_{2}(0)$
$\displaystyle=\frac{S_{2}(0)}{N_{2}},\qquad
e_{2}(0)=\frac{E_{2}(0)}{N_{2}},\qquad i_{2}(0)=\frac{I_{2}(0)}{N_{2}},\qquad
r_{2}(0)=\frac{R_{2}(0)}{N_{2}}.$

##  3 Parameter Identification

In the following, we outline the methodology for calibrating the model
parameters to align closely with observed epidemiological data. This process
is crucial for ensuring that our simulations are able to accurately reflect
the dynamics of a given infectious disease spread. The fitting is based on the
(slightly modified) Levenberg-Marquardt algorithm. We focus on identifying the
parameters of the models using the PDE model, where the number of infectious
individuals serves as our primary target. We specifically analyze the fraction
of population density of infectious individuals $i(x,t;p)$ with parameters
$p=\begin{pmatrix}\sigma&\phi_{e}&\phi_{i}&\beta_{i}\end{pmatrix}^{T}$,
assuming $\beta_{e}=\beta_{i}$ like in [33, 34]. Note that the fraction of
population density of infectious $i$ depends only on three out of four
parameters explicitly, but on all parameters implicitly. Further, let
$\hat{\mathcal{I}}_{k}^{(l)}$ be the target value of the number of infectious
individuals

$\displaystyle\mathcal{I}^{(l)}(t_{k})=\int_{\Omega^{(l)}}i(x,t_{k})n(x)\;dx\in\mathbb{R}$

at time $t_{k}$, $k=1,\dots,K$ for subdomain $\Omega^{(l)}$ with
$l=1,\dots,L$, where $i(x,t_{k})$ is the fraction of population density of
infectious at location $x$ at time $t_{k}$ with $k=1,\dots,K$.

Applying Levenberg-Marquardt algorithm (see Algorithm 1 in Appendix A, similar
to [15, 22, 37]), we have to solve

$\displaystyle\bigl{(}\text{DF}(p^{(l)})^{T}\text{DF}(p^{(l)})+\lambda_{1}I\bigr{)}\Delta
p^{(l)}=-\text{DF}(p^{(l)})^{T}F(p^{(l)})$

with identity matrix $I$, first damping coefficient $\lambda_{1}\geq 0$, new
parameters

$\displaystyle p^{(l+1)}$ $\displaystyle=p^{(l)}+\lambda_{2}\Delta p^{(l)}$

with second damping coefficient $\lambda_{2}\in(0,1]$, error function

$\displaystyle F(p)$
$\displaystyle=\bigl{(}\mathcal{I}^{(l)}(t_{k};p)-\hat{\mathcal{I}}_{k}^{(l)}\bigr{)}_{k=1,l=1}^{K,L}\in\mathbb{R}^{KL}$

and their Jacobian matrix

$\displaystyle\text{DF}(p)$
$\displaystyle=\biggl{(}\frac{\partial\mathcal{I}^{(l)}(t_{k};p_{m})}{\partial
p_{m}}\biggr{)}_{k=1,l=1,m=1}^{K,L,M}\in\mathbb{R}^{KL\times M},$

where $p=\bigl{(}p_{m}\bigr{)}_{m=1}^{M}$ are our $M=4$ parameters we want to
optimize. The damping coefficient $\lambda_{1}$ is chosen so that the matrix
$\bigl{(}\text{DF}(p^{(l)})^{T}\text{DF}(p^{(l)})+\lambda_{1}I\bigr{)}$ is
strictly diagonally dominant. Since the matrix
$\bigl{(}\text{DF}(p^{(l)})^{T}\text{DF}(p^{(l)})+\lambda_{1}I\bigr{)}$ is
symmetric and has real values, we therefore get a positive definite matrix.
For $\lambda_{1}=0$ we apply the Gauss-Newton method and for $\lambda_{1}\gg
0$ the gradient descent method [22].

Now, we want to determine
$\frac{\partial\mathcal{I}^{(l)}(t_{k};p_{m})}{\partial p_{m}}$. We know that

$\displaystyle\frac{\partial\mathcal{I}^{(l)}(t_{k};p_{m})}{\partial p_{m}}$
$\displaystyle=\int_{\Omega^{(l)}}\frac{\partial\bigl{(}i(x,t_{k};p_{m})n(x;p_{m})\bigr{)}}{\partial
p_{m}}\;dx$ $\displaystyle=\int_{\Omega^{(l)}}n(x;p_{m})\frac{\partial
i(x,t_{k};p_{m})}{\partial p_{m}}+i(x,t_{k};p_{m})\frac{\partial
n(x;p_{m})}{\partial p_{m}}\;dx$
$\displaystyle=\int_{\Omega^{(l)}}n(x;p_{m})\frac{\partial
i(x,t_{k};p_{m})}{\partial p_{m}}\;dx,$

i.e., we need to determine the derivative of the fractional population density
$i$ with respect to the parameters $p_{m}\;\forall m=1,\dots,M$. Recall our
PDE for infectious

$\displaystyle\partial_{t}i$
$\displaystyle=\mathop{\mathrm{missing}}{div}\bigl{(}D\nabla i\bigr{)}+\sigma
e-\phi_{i}i$ $\displaystyle=\mathop{\mathrm{missing}}{div}\bigl{(}D\nabla
i\bigr{)}+p_{1}e-p_{3}i.$

Let $f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)\coloneqq\mathop{\mathrm{missing}}{div}\bigl{(}D\nabla
i\bigr{)}+p_{1}e-p_{3}i$ be the right hand side of the PDE for infectious with
the solution of the PDE model
$y=\bigl{(}y_{n}\bigr{)}_{n=1}^{N}=\begin{pmatrix}s&e&i\end{pmatrix}^{T}\in\mathbb{R}^{3}$.
Differentiating $f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)$ with respect to $p_{m}$, $m\in\\{1,\dots,4\\}$ results
in

$\displaystyle\partial_{p_{m}}\bigl{(}\partial_{t}i\bigr{)}$
$\displaystyle=\partial_{p_{m}}\bigl{(}f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)\bigr{)}$
$\displaystyle=\sum_{n=1}^{N}\partial_{\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})}f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)\cdot\partial_{p_{m}}\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})\bigr{)}$
$\displaystyle\qquad+\partial_{y_{n}}f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)\cdot\partial_{p_{m}}y_{n}+\partial_{p_{m}}f_{i}(\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})_{n=1}^{N},y;p)\in\mathbb{R}$

for $m=1,\dots,M$ with

$\displaystyle\partial_{p_{m}}\bigl{(}\partial_{t}i\bigr{)}$
$\displaystyle=\sum_{n=1}^{N}\partial_{\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})}f_{i}\cdot\partial_{p_{m}}\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})\bigr{)}+\partial_{y_{n}}f_{i}\cdot\partial_{p_{m}}y_{n}+\partial_{p_{m}}f_{i}$
$\displaystyle=\underbrace{\partial_{\mathop{\mathrm{missing}}{div}(D\nabla
i)}f_{i}}_{=1}\cdot\partial_{p_{m}}\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
i)\bigr{)}+\underbrace{\partial_{s}f_{i}}_{=0}\cdot\partial_{p_{m}}s+\underbrace{\partial_{e}f_{i}}_{=p_{1}}\cdot\partial_{p_{m}}e+\underbrace{\partial_{i}f_{i}}_{=-p_{3}}\cdot\partial_{p_{m}}i+\partial_{p_{m}}f_{i}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla\bigl{(}\partial_{p_{m}}i\bigr{)})+p_{1}\partial_{p_{m}}e-p_{3}\partial_{p_{m}}i+\delta_{m,1}e-\delta_{m,3}i$

in terms of Kronecker’s delta $\delta_{j,k}$.

Now, considering the differential equation for exposed

$\displaystyle\partial_{t}e=f_{e}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
e)+\biggl{(}1-\frac{A}{n_{l}+n_{0}}\biggr{)}s(\beta_{i}i+\beta_{e}e)-\sigma
e-\phi_{e}e$ $\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
e)+\biggl{(}1-\frac{A}{n_{l}+n_{0}}\biggr{)}p_{4}s(i+e)-p_{1}e-p_{2}e,$

we have to solve the following DEs

$\displaystyle\partial_{p_{m}}\bigl{(}\partial_{t}e\bigr{)}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla\partial_{p_{m}}e)+\left(1-\frac{A}{n_{l}+n_{0}}\right)\big{(}\delta_{m,4}s(e+i)+p_{4}(e+i)\partial_{p_{m}}s+p_{4}s(\partial_{p_{m}}e+\partial_{p_{m}}i)\big{)}$
$\displaystyle\quad-
p_{1}\partial_{p_{m}}e-p_{2}\partial_{p_{m}}e-\delta_{m,1}e-\delta_{m,2}e.$

Moving on to the DE for susceptibles

$\displaystyle\partial_{t}s=f_{s}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
s)-\biggl{(}1-\frac{A}{n_{l}+n_{0}}\biggr{)}s(\beta_{i}i+\beta_{e}e)$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla
s)-\biggl{(}1-\frac{A}{n_{l}+n_{0}}\biggr{)}p_{4}s(i+e),$

we must solve the subsequent DE

$\displaystyle\partial_{p_{m}}\bigl{(}\partial_{t}s\bigr{)}$
$\displaystyle=\mathop{\mathrm{missing}}{div}(D\nabla\partial_{p_{m}}s)-\bigl{(}\mathop{\mathrm{missing}}{div}(D\nabla
y_{n})\bigr{)}\big{(}\delta_{m,4}s(e+i)+p_{4}(e+i)\partial_{p_{m}}s+p_{4}s(\partial_{p_{m}}e+\partial_{p_{m}}i)\big{)}.$

As the initial conditions of the PDE system (1) do not depend on the
parameters $p$, the initial conditions of the variational equations [10] or
sensitivity equations [34] with solution
$\begin{pmatrix}\partial_{p_{m}}s&\partial_{p_{m}}e&\partial_{p_{m}}i\end{pmatrix}^{T}$
are zero. The boundary conditions of the sensitivity equations are the same as
for the PDE system (1). We solve the system of PDEs and variational equations
simultaneously. This procedure is often termed as forward sensitivity analysis
[24].

##  4 Implementation

In this section, we detail the technical processes and computational methods,
and tools underlying the implementation of our hybrid epidemiological model.

To start, ensuring positive initial conditions for the PDE system is
foundational. We impose the constraint $w_{y_{n}}\geq 0$ for $y_{n}=s,e,i,r$
and solve the corresponding constraint optimization problem using Sequential
Least Squares Programming (SLSQP, as implemented in SciPy 1.12.0 using the
scipy.optimize.minimize function).

Building on this initial setting, we proceed to calibrate the initial
parameters, which is crucial for capturing the nuanced dynamics of disease
spread. This calibration is achieved by fitting the fractional infectious
density of the derived ODE model by using different patches, each patch
representing one subdomain. Notably, inter-patch transmission is excluded from
this phase to simplify the fitting process. This fitting step is done by
minimizing the norm between target values and our numerical solution. The
numerical solution is computed using finite differences with time step size
$\Delta t=0.1$, again by using SLSQP.

Further refining the model, parameter identification for the PDE components
employs a slightly modified version of the Levenberg-Marquardt algorithm, as
detailed in Appendix A. This step aims to minimize the weighted error
$e_{w}^{(l)}=||w_{F}(\Omega)F(p^{(l)})||_{2}$ where the weights vary by domain
importance, amplifying the significance of Very Important Domains (VIDs), such
as Milan (see below). The weight vector $w_{F}(\Omega)$ is defined as:

$\displaystyle w_{F}(\Omega)=\begin{cases}2,&\text{if }\Omega\text{ is
VID,}\\\ 1,&\text{if }\Omega\text{ is not VID.}\end{cases}$

This weighted error strategy allows for simultaneous optimization of all
parameters and also a focused adjustment of specific parameters, such as
$\beta_{e}$. This approach is inspired by [33, 34]. The optimization continues
iteratively until the error falls below a predefined threshold or the maximum
number of iterations is reached. Then, we proceed with the subsequent temporal
segments.

To support these computational efforts, the triangulation of the modeling
domain was performed using Triangle [29], while the full-PDE and the hybrid
PDE-ODE systems were solved using the finite element method implemented in the
Kaskade7 software [18], using the Dune interface. Finally, the visualization
of our results was conducted using ParaView [31], a tool well-suited for
rendering and analyzing complex simulation data.

##  5 Numerical Results

In this section, we focus on a detailed exploration of our hybrid model’s
efficacy through a series of numerical experiments. Our aim is to evaluate the
performance of the hybrid model in comparison to the conventional full-PDE
model across a spectrum of scenarios, ranging from controlled synthetic
environments to complex real-world settings. To lay the groundwork, we begin
by selecting a simplistic geographical area as our testing ground, where we
vary the division of the domain into PDE and ODE segments (see section 5.1).
Here, we aim not just to compare the hybrid model with its full-PDE
counterpart, but also to illustrate the trade-offs between computational
efficiency and solution accuracy (see section 5.1.1). This initial phase is
crucial for quantifying the model’s stability and its capacity to maintain
precision across different configurations, thereby providing us with a
foundational understanding of the model’s capabilities.

Moving beyond synthetic experiments, we advance to the evaluation of the
models against empirical data, bringing into focus the real-world
applicability and validity of our modeling approach. Our investigation extends
across two distinct locales: Lombardy, Italy, and Berlin, Germany for which
real-world data is available (see section 5.2). This validation approach not
only tests the models’ predictive precision but also their adaptability and
responsiveness to real-world infectious disease dynamics. The results of these
experiments will help us understand how our hybrid modeling approach can be
used in the complex field of epidemiological forecasting and analysis.

###  5.1 Experiments Based on Synthetic Data

In these set of experiments, we investigate the relationship between solution
accuracy and computational efficiency within a partitioned rectangular domain
$\Omega=\mathopen{]}0,2\mathclose{[}\times\mathopen{]}0,1\mathclose{[}$,
dividing it into PDE and ODE regions with varying percentages. Our aim is to
uncover potential trade-offs in achieving accuracy while maintaining
computational efficiency. Additionally, we explore extreme cases within a
50-50 domain division, assessing the model’s resilience when infected
individuals are concentrated in one region and evaluating its response to
boundary challenges. These experiments provide valuable insights into the
model’s behavior under different configurations, contributing to a
comprehensive understanding of the system’s performance.

####  5.1.1 Trade-offs Between Solution Accuracy and Computational Efficiency

The PDE model’s solution is assumed to represent an idealized benchmark
against which the results from the ODE model and the hybrid PDE-ODE model are
compared. Considering the entire domain for evaluation, any error arising from
the PDE domain or the region $\Omega_{1}\subset\mathbb{R}^{2}$ of interest,
modeled in spatial detail, is already incorporated through the coupling. The
error measurement is conducted as an absolute error, calculated as the sum of
deviations from the benchmark. The accuracy is determined using the formula

$\displaystyle\text{accuracy}=1-\frac{||\text{deviation}||_{2}}{\text{maximum
possible deviation}}.$

Here, the solution of the PDE model achieves an accuracy of 1, while an ODE
model’s solution attains an accuracy of 0, effectively reflecting the maximum
possible deviation.

percentage of ODE-domain | number of
---|---
| vertices | edges | triangles
0$0$ | $25,669$ | $71,004$ | $45,336$
0$12.5$ | $22,918$ | $63,251$ | $40,334$
$25$ | $20,113$ | $55,336$ | $35,224$
0$37.5$ | $17,336$ | $47,505$ | $30,170$
$50$ | $14,561$ | $39,680$ | $25,120$
0$62.5$ | $11,753$ | $31,756$ | $20,004$
$75$ | 0$8,941$ | $23,820$ | $14,880$
0$87.5$ | 0$6,177$ | $16,028$ | 0$9,852$
$100$ 0 |  $0$ 00 |  $0$ 00 |  $0$ 00
Table 2: Information on mesh of PDE-domain of the hybrid model (4) and full-
PDE model (1).

We want to partition the domain vertically and vary the sizes for the PDE-
domain and the ODE-domain. We simultaneously adjust the proportion of infected
and recovered individuals in their respective regions, along with the total
population, during the proportional partitioning of the area. We scale the
infected individuals to represent the exposed category. The fractional
densities $s$, $e$, $i$ are initially set as constants across the region, and
for the population density, we use Gaussian functions that remain constant
along the x-axis with a variance of $0.1$. With this choice of the Gaussian
function, we prevent disadvantages that may arise when partitioning the
domain. In Table 3, you will find the parameters specific to the rectangular
domain across various time intervals.

parameter | value
---|---
| time index 0 | time index 10 | time index 23 | time index 29
$\sigma$ | $2.6676e-02$ | - | - | $1.8747e-01$
$\phi_{e}$ | $0.0000e+00$ | - | - | $3.2655e-11$
$\phi_{i}$ | $2.3310e-01$ | - | - | $2.1907e+00$
$\beta_{i}$, $\beta_{e}$ | $4.4202e-01$ | $2.0588e-01$ | $6.0352e-08$ | $2.2168e-01$
$D$ | $\tfrac{0.0435}{111.3^{2}}\cdot 1.0e+02$ | $\tfrac{0.0198}{111.3^{2}}\cdot 1.0e+2$ | $\tfrac{0.0090}{111.3^{2}}\cdot 1.0e+2$ | $\tfrac{0.0075}{111.3^{2}}\cdot 1.0e+2$
$A$ | $3.0e+07$ | - | - | -
$n_{0}$ | $4.5e+07$ | - | - | -
Table 3: Parameters of PDE model (1) and of the hybrid model (4) for
rectangular domain. The symbol “-” means that the parameter has not changed
over time.

The infectious density over time can be observed in the following Fig. 1.
Given that the fractional densities $s$, $e$, $i$ are initially constant and
the population density $n$ follows a Gaussian distribution, the resulting
infectious density $i\cdot n$ will also follow a Gaussian distribution for
$t=0$. The Gaussian-shaped spread of infectious disease is evident due to the
constant Gaussian-shaped population density.

Figure 1: Infectious density of PDE model (1) for times $t\in\\{0,29,30,59\\}$
in rectangular domain. The Gaussian-shaped spread of infectious disease is
evident due to the Gaussian-shaped population density.

Figure 2 shows the number of infectious individuals for each time step, along
with the maximum average deviation in the number of infectious individuals
between the hybrid model and the full-ODE model when compared to the full-PDE
model. The maximum average deviation measures imply that the variations
between hybrid models are around 0.99%, as calculated by the differences in
the right plot. In other words, a 12.5% difference in area corresponds to an
approximate 1% error when comparing hybrid models. A more notable error
difference is observed for the hybrid model in the right plot with an ODE
domain size of 87.5% compared to the full-ODE model, which is around
$7.95\%-7.21\%=0.74\%$. Overall, these error differences suggest that the
hybrid models tend to shift towards the direction of the ODE model, indicating
that the ODE component of the hybrid model seems to exert a stronger weighting
or influence on the boundary conditions between the regions, consequently
affecting the model results of the PDE component of the hybrid model to a
greater extent.

Figure 2: Number of infectious individuals of full-PDE and hybrid model in
rectangular domain for times $t\in[0,59]$ (left) and maximum average deviation
from full-PDE model (right). The infectious count decreases with a higher
percentage of the ODE domain. This is attributed to the inclusion of the Allee
term, which spatially modifies the infection rate in the PDE domain while
remaining constant in the ODE domain.

Measuring the absolute deviation to the full-PDE model (see Fig. 3 and Table
4), we can observe an approximately linear growth. If the focus is on
enhancing predictions or attaining greater precision in the results, it is
recommended to choose the full-PDE model or consider a more advanced model,
which may include additional compartments [14] or even an ABM [16]. We assume
that the runtime scales proportionally with the PDE domain size, considering
the ODE system can be solved in a fraction of the time and is regarded as a
non-critical contributor to the duration. In this context, when comparing the
program runtime with the accuracy of the results, it becomes evident that
sacrificing a small amount of accuracy in favor of solving the model more
efficiently is a viable option.

It is noticeable that the infectious count decreases with a higher percentage
of the ODE domain, a trend typically not observed [4], [28]. This is
attributed to the inclusion of the Allee term, which spatially modifies the
infection rate in the PDE domain while remaining constant in the ODE domain.
Setting $A=0$ results in almost identical outcomes for all models, which is
why we deliberately retained the Allee term.

Figure 3: Accuracy (mean absolute error) of the full-PDE model (1) and hybrid
model (4) in rectangular domain. We can observe an approximately linear
growth.
percentage of ODE-domain | error
---|---
0$0$ | 0$0.00$
0.$12.5$ | 0$3.30$
$25$ | 0$5.95$
0.$37.5$ | 0$8.61$
$50$ | $11.27$
0.$62.5$ | $13.93$
$75$ | $16.59$
0.$87.5$ | $19.25$
$100$ | $21.27$
Table 4: Accuracy (mean absolute error) of the full-PDE model (1) and hybrid
model (4).

####  5.1.2 Extreme Cases

Exploring extreme cases within a 50-50 domain division, we assess the model’s
resilience by adjusting the total number of infectious individuals in the PDE
contribution of the hybrid model to match the total population number.
Simultaneously, we set it to zero in the ODE, and vice versa. Additionally, we
introduce a Gaussian circular function that remains constant along the y-axis
directly at the boundary $\Gamma$ at $x=1.0$. The results are shown in Figure
4. The results for each domain are shown in Figure 4. Due to significant
differences between the models, the infectious densities of the two models
cannot be effectively visualized for comparison. Simulation plots of the
fractional infectious density for the hybrid model were created for an extreme
case (see Fig. 5), serving the purpose of visualization. These simulation
plots provide no additional value, hence, the exclusion of spatial
visualization for other extreme cases. It is evident that the outcomes for the
ODE region appear visually unaffected by the location of the population in the
PDE domain. However, it is noticeable that there is a flow of infections into
the ODE, indicating that the boundary conditions influence the ODE. This is
evident as infection numbers increase in the right plot from zero to values as
high as 8,256.5 (yellow) and 8,246.43 (green). Due to their simplicity, ODE
models may exhibit a higher sensitivity to specific changes in parameters
compared to PDE models. Starting at time $t=10$, a new time interval begins,
introducing updated parameters. Visually, this adjustment has no apparent
impact on the infection dynamics within the PDE, in contrast to the ODE, which
experiences a second minor wave. This wave rapidly diminishes with the
adoption of the new parameters at time $t=29$, leading to a noticeable decline
in infection numbers within the PDE region as well.

Figure 4: Extreme cases of the hybrid model (4) in rectangular domain: the
number of infectious individuals is initially equal to the total population
number (left) and zero (right). The outcomes for the ODE region appear
visually unaffected by the location of the population in the PDE domain,
whereas the outcomes for the PDE region vary significantly.

Figure 5: Fractional infectious density of hybrid model (4) for times
$t\in\\{0,2,10,29\\}$ in rectangular domain. Initially, the number of
infectious individuals for the PDE contribution of the hybrid model is equal
to zero, while for the ODE contribution, it is equal to the population number.
Due to the boundary conditions between the two regions, there is an infection
flow from the ODE region (right) to the PDE region (left).

When examining the right plot of Figure 4, where the initial number of
infectious individuals is zero, the PDE with horizontal population density
(red) exhibits a significantly lower peak in the infection curve compared to
the PDE with vertical population density (blue). Due to the boundary
conditions between the two regions, there is an infection flow from the ODE
region to the PDE region (see e.g. Fig. 5). This is more pronounced because
the population density is increased precisely at the boundary and diminishes,
even approaching zero, in the rest of the region. Consequently, the infection
flow has a stronger impact. In the case of horizontal population density,
either the population density is very low where the fractional infectious
density is higher, or the fractional infectious density is low where the
population density is high. As a result, the infection flow has a lesser
effect.

A similar principle applies to the PDE results in the left plot. When
examining the left plot (and the corresponding data), where the initial number
of infectious individuals corresponds to the population size, the PDE with
horizontal population density (green) has fewer total infectious than the ODE
(red, blue) for all time steps except the last few. Due to the boundary
conditions and the absence of infectious in the ODE, the fractional infectious
density of the PDE is reduced at the boundary, once again exerting a stronger
influence on the vertical population density than the horizontal. Meanwhile,
the infection numbers in the ODE region increase, and this, in turn, has a
more pronounced impact on the rising infection numbers in the PDE for the
vertical population density (yellow) than the horizontal (green).

Understanding the infection dynamics in a hybrid PDE-ODE model with Dirichlet
boundary conditions is crucial for unraveling the interplay between different
regions. The chosen boundary conditions significantly influence the infection
dynamics within this hybrid framework. For instance, if the fractional
infectious density in the ODE region is higher than at the boundary of the PDE
region, the simulation exhibits a migration of infected individuals from the
interior of the PDE region toward the boundary to fulfill the Dirichlet
conditions. This leads to an accumulation of infected individuals at the
boundary. Conversely, in cases where the fractional infectious density in the
ODE region is lower than at the boundary of the PDE region, the system
demonstrates a diffusion of infected individuals into the interior of the PDE
region. In such scenarios, the infection dynamics would be predominantly
influenced by the PDE model, with minimal impact observed in the ODE model.

To amplify the influence of the Dirichlet boundary conditions and observe
heightened effects in the ODE region, adjusting the penalty term could be
employed. By assigning a higher weight to the boundary conditions through an
elevated penalty term, the model emphasizes the importance of meeting the
specified conditions at the boundary, potentially leading to more pronounced
effects in the ODE solution.

###  5.2 Real-World Examples

In the following section, we evaluate the practical applicability of our
hybrid PDE-ODE model through its application to actual epidemiological
scenarios. This examination not only validates the model against empirical
data but also highlights its potential for real-world utility in public health
management. We illustrate the model’s effectiveness by applying it to two
distinct geographical settings: Lombardy, Italy, and Berlin, Germany. These
locations were chosen due to their varied population densities and distinct
public health challenges during the COVID-19 pandemic, offering a robust
framework for testing and demonstrating the model’s capabilities. For each
region, we carefully calibrated the model using available epidemiological
data, adapting our simulations to reflect the specific disease dynamics
observed in each area. The results from these case studies provide insights
into the model’s ability to simulate and predict the complexities of
infectious disease dynamics in the context of diffusive processes.

####  5.2.1 Lombardy, Italy

In the first example, we explore the application of our hybrid PDE-ODE model
to the Italian region of Lombardy (see Fig. 6). We follow the approach by
Viguerie et. al [33] and focus on the five largest provinces, which were
mainly affected by the pandemic. We compare the performance of our hybrid
model with a full-PDE model. Specifically, in the hybrid model, the province
of Milan, which is relatively small but has a high population, is simulated
using an ODE due to its high population density, facilitating faster and more
uniform population exchanges, making it a more suitable candidate for ODE
simulation compared to the other provinces. The aim of this example is to
demonstrate the predictive power of our model with respect to the spread of
Covid-19. The geographical data, i.e. (1) latitude and longitude of provinces,
(2) boundary points of Lombardy, and (3) boundary points of Milan were taken
from public repositories [1, 2]111Note that in the case of longitude, Lombardy
is close to the equator, i.e., we can simply scale the location data by
$111.3\;\nicefrac{{\text{degrees}}}{{\text{km}}}$ to convert [3]. The same
holds for the latitude, since each degree of latitude is approximately $111.3$
km apart. This is not necessary in the calculations and will ultimately occur
before the dynamics are simulated.. We calibrated our model using documented
case data from Lombardy [1], calculating a 7-day average from these cases as
target for parameter identification.

Figure 6: Main provinces affected by the pandemic in Lombardy: Lodi, Milan,
Cremona, Bergamo, Brescia.

To simulate the dynamics within this region, the mesh setup consists of 63,017
triangles, 99,084 edges and 36,068 vertices. When adapted for the coupled
model, this results in a mesh configuration of 31,953 triangles, 52,912 edges,
and 20,960 vertices. For the initial conditions, we use circular Gaussian
functions with a variance of $0.05$. The Lombardy specific model parameters
are listed in Table 5 and are adjusted across four distinct time intervals.
Start- and end-times are taken from Viguerie et al. [33]. Each interval
corresponds to different phases of pandemic response, reflecting changes in
public health policies and restrictions. These parameters are finely tuned for
each period and remain constant over the time interval. During the initial
time interval, no restrictions are in effect. However, as each subsequent time
interval begins, restrictions are introduced, and changes in parameter values
show how policy changes affect disease dynamics.

parameter | value
---|---
| February 27th, 2020 | March 9th, 2020 | March 22nd, 2020 | March 28th, 2020
$\sigma$ | $2.6676e-02$ | - | - | $1.8747e-01$
$\phi_{e}$ | $0.0000e+00$ | - | - | $3.2655e-11$
$\phi_{i}$ | $2.3310e-01$ | - | - | $2.1907e+00$
$\beta_{i}$, $\beta_{e}$ | $4.4202e-01$ | $2.0588e-01$ | $6.0352e-08$ | $2.2168e-01$
$D$ | $\tfrac{0.0435}{111.3^{2}}\cdot 1.0e+02$ | $\tfrac{0.0198}{111.3^{2}}\cdot 1.0e+2$ | $\tfrac{0.0090}{111.3^{2}}\cdot 1.0e+2$ | $\tfrac{0.0075}{111.3^{2}}\cdot 1.0e+2$
$A$ | $3.0e+07$ | - | - | -
$n_{0}$ | $4.5e+07$ | - | - | -
Table 5: Parameters of PDE model (1) and of the hybrid model (4) for Lombardy.
The symbol “-” means that the parameter has not changed over time.

Figure 7: Fractional infectious density in km2 of PDE model (1) (left) and PDE
contribution of hybrid model (4) (right) in Lombardy for times
$t\in\\{1,28,59\\}$. The boundary conditions between the PDE and ODE domains
ensure, among other things, a uniform flow of infectious originating from
Milan in all directions.

Figure 8: Number of infectious people of PDE (1) and hybrid model (4) in each
province (simulated and real 7-day average data). In Cremona, Bergamo, and
Brescia, both models showed similar trends, but the hybrid model overestimated
infections in Milan and underestimated them in Lodi.

The dynamics are simulated and shown in Fig. 7, and the total number of
infectious individuals are depicted in Fig. 8. A notable feature of our
simulations is the incorporation of inter-provincial interactions,
particularly the bidirectional transmission of infections between Lodi and
Milan, facilitated by specified boundary conditions. This dynamic interaction
results in a significant impact on infection trends in both provinces.

Our findings revealed a consistent pattern across both modeling approaches in
provinces Cremona, Bergamo, and Brescia (see Fig. 8). The simulation results
for these provinces showed similar trends, indicating that both models provide
consistent and reliable predictions, regardless of the chosen modeling
approach, for these provinces located farther from the boundary of the ODE-
modeled area. However, a noteworthy deviation emerged in the hybrid model’s
predictions for Milan and its neighboring province, Lodi (see Fig. 8). The
hybrid model tended to overestimate the infectious cases in Milan.
Simultaneously, it underestimated the infection rates in the adjacent province
of Lodi. The boundary conditions between the PDE and ODE domains ensure, among
other things, a uniform flow of infectious originating from Milan in all
directions, a feature that is clearly visible in Fig. 7. This leads to a
slight increase in infectious numbers in Bergamo. Simultaneously, the
infectious flow from Lodi towards Milan is restrained due to the boundary
conditions. Instead, infectious cases tend to migrate more prominently towards
the neighboring province of Cremona.

The parameter fitting for the full-PDE model naturally yields results closer
to the target data than those obtained from the hybrid model. An alternative
approach involves directly fitting the parameters of the hybrid model.

h

h

####  5.2.2 Berlin, Germany

For the second example, we apply our hybrid PDE-ODE model to the urban context
of Berlin, Germany. Berlin presents an interesting case with its diverse and
dense urban structure, which offers a challenging environment for
epidemiological modeling due to varied population distribution and mobility
patterns.

The geographical data, i.e. latitude and longitude of Berlin, is publicly
available222http://opendatalab.de/projects/geojson-utilities/. Utilizing this
data, as before, we created a mesh for the full-PDE model, consisting of
73,294 triangles, 110,471 edges and 37,178 vertices. For our hybrid model, we
divided the domain vertically. This resulting reduced domain consists of
36,633 triangles, 55,249 edges and 18,617 vertices.

However, opposed to the Lombardy study, this time we used data generated by an
agent-based model (ABM) specifically calibrated for Berlin (see [25] for
details). This ABM effectively simulates individual-level behaviors that
influence disease transmission dynamics. The model is based entirely on
mobility data, derived from mobile phone data, and incorporates a 25%
sample333https://doi.org/10.14279/depositonce-11495.2 of the Berlin
population. This allows to capture typical movements of persons as well as
their activity reductions during the various public health restrictions. While
the foundational ABM parameters were sourced from existing literature, the
infection rate parameter was specifically calibrated to align with the
observed data during the COVID-19 pandemic.

While Lombardy’s model focused on the dynamics across larger provincial areas,
the Berlin model incorporates the granularity of individual agents,
representing residents with specific demographic and health status
characteristics. This approach allows us to capture detailed spatial
information while maintaining a constant population density over time,
providing a comprehensive understanding of infection dynamics in urban
environments, such as Berlin.

Each agent is assigned home coordinates (see Fig. 9). We map these home
coordinates to grid points based on the shortest distance. This enables us to
determine the number of agents at each grid point with specific health
statuses and how many agents are present there. Following this, we have
smoothed the (sub)populations on the mesh. Subsequently, we calculate the
density at each grid point. Initially, we determine the population density
$n$, followed by the fractional densities $s$, $e$, $i$, $r$. Since not all
locations have residences in our area, resulting in zero population density in
those areas, we fill the gaps with the smallest representable value to
minimize their impact on the overall population. Filling the gaps ensures our
condition $s+e+i+r=1$ holds. Additionally, this population fraction is labeled
as removed to ensure that this additional population does not affect the
infection dynamics.

Figure 9: Original home coordinates of agents in Berlin, Germany (left) and
adjusted locations to grid (right).

Although Berlin is not situated close to the equator, for the purpose of our
simulation, we have chosen a simplified approach by scaling location data
using $111.3\;\nicefrac{{\text{degrees}}}{{\text{km}}}$. While this method may
introduce some margin of error due to the deviation from the equator, it
allows for a straightforward representation of spatial relationships in our
model. The restrictions were introduced in Berlin at different times and to a
different extent than in Lombardy [8]. We chose March 15th, 2020, and March
23rd, 2020444https://taz.de/Corona-Tagebuch-der-Hauptstadt/!5671055/, as the
pivotal dates for the introduction of (additional) restrictions. We could not
find details about the last restriction in Lombardy and, therefore, omitted it
in the case of Berlin. In Table 6, you can find the parameters corresponding
to different time intervals for Berlin.

parameter | value
---|---
| February 25th, 2020 | March 15th, 2020 | March 23rd, 2020
$\sigma$ | $1.8513e+01$ | $1.8513e+01$ | $1.8513e+01$
$\phi_{e}$ | $2.9179e-03$ | $1.8561e-02$ | $8.3704e-03$
$\phi_{i}$ | $2.7734e-01$ | $3.1585e-01$ | $2.8986e-01$
$\beta_{i}$, $\beta_{e}$ | $2.9740e+00$ | $1.4135e+00$ | $1.4490e+00$
$D$ | $\tfrac{0.0435}{111.3^{2}}\cdot 1.0e+02$ | $\tfrac{0.0198}{111.3^{2}}\cdot 1.0e+2$ | $\tfrac{0.0090}{111.3^{2}}\cdot 1.0e+2$
$A$ | $1.0e+09$ | - | -
$n_{0}$ | $1.5e+09$ | - | -
Table 6: Parameters of PDE model (1) and of the hybrid model (4) for Berlin.
The symbol “-” means that the parameter has not changed over time.

Figure 10: Infectious density of PDE model (1) for times
$t\in\\{1,8,18,60\\}$. Here, the spread can be observed particularly well at
the level of home locations.

We adopted the diffusion coefficient for Berlin in a manner consistent with
that used for Lombardy, although it’s worth noting that different time
intervals were employed, and Berlin was simulated for an overall longer
duration. However, due to Berlin’s significantly smaller area, the spread in
Berlin appears faster in relation to the geographical size. This is clearly
visible in Figs. 11, 14, where we have plotted the fractional infectious
density. The spread of infections can be observed particularly well at the
level of home locations in Fig. 10. Unsurprisingly, once again we observe that
the fitting of the PDE model contributes to the closer alignment with the
target data (see Fig. 12). We observe a quite similar issue here, which we
have already seen in experiments conducted in a rectangular area (see Sec.
5.1.2). The uniform flow of infection from the ODE region, combined with a
high population density directly at the boundary between the differently
modeled regions, dominates, even though the population density is also high in
other areas of Berlin (compare to susceptible density in Fig. 16). This
results in a sharp increase in infection numbers in the hybrid model. These
two exposed agents in the north have a stronger effect on the spread than the
other two exposed agents further south (see Fig. 15). This is likely because
we have no-flux Neumann boundary conditions in the north, preventing the
infection from spreading out of Berlin. Note that the population density
roughly corresponds to the susceptible density initially, and in Figure 16, it
can be seen that the population density for the exposed agents in the north is
lower than those further south. Directly at the boundary, the fractional
infectious density is slightly higher than in the adjacent regions. This is
again due to the Dirichlet boundary conditions we enforce adjacent to the ODE
region. The infection flow slowly reaches the densely populated region in
western Berlin, as we have a large area in the region with very low population
density, and many areas between the home locations with very low population
density, thus resulting in lower infection rates in these areas. Additionally,
we only have an infection flow originating from the boundary between the PDE
and ODE regions, along with a few initial COVID cases that are not located in
western Berlin. Including infection flows from outside Berlin would be
advantageous here. Once the flow across the boundary is prohibited, by setting
the boundary conditions between the two regions to zero, it becomes apparent
that the infectious numbers only explode within the ODE region (see Fig. 17).
Originally, the boundary conditions were intended to maintain the fractional
densities on both sides of the boundary as similar as possible, which is why
the fractional infectious density in the PDE region has increased
significantly. This, in turn, can cause an escalating effect in the ODE
region, as we can observe. The primary reason for the varying number of
infectious cases is likely due to the non-uniform distribution of the
population within the PDE domain, with gaps between the home locations.
Achieving a uniform distribution of home locations is not feasible in this
experiment, as we specifically aimed to test precise data and its impact on
the infection dynamics of the hybrid model. What we did in this regard was to
assign home locations to a grid point and distribute them more evenly within
the adjacent triangles whose vertices touch the very same grid point.
Additionally, one could consider a non-adaptive coarser grid, which would make
the solution less accurate. Further experiments are required in this regard
and will be explored in the future.

Figure 11: Fractional infectious density of PDE model (1) for times
$t\in\\{1,8,18,60\\}$. Due to Berlin’s significantly smaller area, the spread
in Berlin appears faster than in Lombardy.

Figure 12: Number of infectious people of full-PDE model (1) (left) and hybrid
model (4) (right) in Berlin (simulated and ABM data). The primary reason for
the varying number of infectious cases is likely the non-uniform distribution
of the population within the PDE domain.

Figure 13: Infectious density of hybrid model (4) for times
$t\in\\{1,8,18,60\\}$. Here, the spread can be observed particularly well at
the level of home locations.

Figure 14: Fractional infectious density of hybrid model (4) for times
$t\in\\{1,8,18,60\\}$. Due to Berlin’s significantly smaller area, the spread
in Berlin appears faster than in Lombardy.

Figure 15: (Fractional) Exposed density of hybrid model (4) for times
$t\in\\{0,1,2\\}$. The scale for the initial time step is incorrect, as the
results in the first image exceed the range of the scale. The initial time
step was included solely to illustrate the initial position of exposed
individuals. By time step two, the infection flow originating from the
boundary has already begun to have a greater impact than the exposed
individuals. Figure 16: Susceptible density of hybrid model (4) for time
$t=0$. The population density roughly corresponds to the susceptible density
initially.

Figure 17: Number of infectious people of PDE contribution (left) and ODE
contribution (right) of hybrid model (4) in Berlin (simulated and ABM data)
with zero boundary conditions. Once the flow across the boundary is
prohibited, it becomes apparent that the infectious numbers only explode
within the ODE region.

##  6 Conclusions

In this paper, we introduced a novel hybrid mathematical modeling approach
that couples Partial Differential Equations (PDEs) with Ordinary Differential
Equations (ODEs) to simulate epidemiological processes effectively across
geographic landscapes. Our model exemplifies the integration of the spatially
detailed representation of disease dynamics provided by PDEs with the
computational efficiency of ODEs, allowing for the timely assessment of public
health interventions and potential impacts of infectious diseases across large
populations. Through the implementation and evaluation of this hybrid PDE-ODE
model in a simple setting using synthetic data and in two distinct real-world
settings, specifically Lombardy and Berlin, we explored how interactions
between ODE and PDE regions can significantly influence infection dynamics.
The results from these case studies highlight the model’s adaptability and its
ability to reflect complex epidemiological behaviors in various settings.

Building on these insights, our findings suggest the necessity for strategic
adjustments to minimize potential risks, particularly the steady flow of
infections that has a heightened impact on neighboring areas with high
population density. To manage the spread of infectious diseases more
effectively, we recommend modifying the size and location of ODE regions or
creating buffer zones within PDE domains. For example, it is advisable to
avoid splitting well-connected areas like Berlin at their center. Instead, the
ODE region can be shifted or resized to minimize its impact on densely
populated areas. If direct changes to the ODE regions are not possible or too
complicated, implementing transition zones within the PDE domains could serve
as an effective alternative. These zones can be established without altering
the proportionate sizes of the regions. For instance, in the Lombardy case-
study, a transition zone was indirectly implemented by not initially placing
individuals from Milan in this area. Instead, a few individuals from Lodi or
Bergamo were distributed through the Gaussian distribution in the transition
zone. Furthermore, on the rectangular domain, a more uniform population
distribution led to better control of infection dynamics. However, the
significant rise in infection numbers in the Berlin ODE region highlights the
need for further experiments to refine this strategy.

Directly fitting the parameters of the hybrid model, rather than adapting
those from the PDE model, would enhance its performance by allowing precise
adjustments tailored to the combined dynamics of the ODE and PDE regions. This
adjustment ensures that the infection numbers for both models align closely
with the target data, providing a more accurate reflection of the
epidemiological situation.

The exploration of infection dynamics beyond the confines of Lombardy and
Berlin introduces additional complexities and opportunities for our modeling
approach. By expanding our analysis to include other regions, the hybrid model
proves very useful for simulating extensive geographical areas and
understanding the broader implications of disease spread. Incorporating
dynamics from external areas could further deepen our understanding and
enhance the model’s applicability.

Finally, considering the urban spread dynamics, introducing a location-
dependent diffusion coefficient could control infection flow at the periphery,
potentially speeding up the spread within cities and modifying or replacing
the Allee term’s role. To further enhance our model, we plan to incorporate
transportation dynamics, including trains and airplanes, using stochastic
differential equations with discrete jumps. By expanding the scope of our
hybrid model to include more detailed interactions and broader geographic
applications, we aim to make ongoing contributions to the field of
epidemiological forecasting and analysis, providing more precise and effective
tools for managing public health crises.

## Acknowledgments

We acknowledge Kai Nagel for valuable discussions on the agent-based model. We
also thank Jakob Benjamin Rehmann and Sydney Cornelia Paltra for providing
essential data for the Berlin experiments and for assistance with their
software.

The work on this paper was funded by the German Ministry of research and
education (BMBF) (project ID: 01KX2022A) and by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy via MATH+: The Berlin Mathematics Research Center
(EXC-2046/1, project ID: 390685689).

## References

  * [1] Public github repository “covid-19”.  https://github.com/pcm-dpc/COVID-19.
  * [2] Public github repository “geojson-italy”.  https://github.com/openpolis/geojson-italy/.
  * [3] Navigation System Design, chapter 8, pages 305–354.  John Wiley and Sons, Ltd, 2020.
  * [4] Hamadjam Abboubakar, Reinhard Racke, and Nicolas Schlosser.  A reaction-diffusion model for the transmission dynamics of the coronavirus pandemic with reinfection and vaccination process.  Technical Report 409, 2023.
  * [5] Carlos J.S. Alves, Nuno F.M. Martins, and Svilen S. Valtchev.  Domain decomposition methods with fundamental solutions for helmholtz problems with discontinuous source terms.  Computers and Mathematics with Applications, 88:16–32, 2021.  Trefftz Methods and Method of Fundamental Solutions - theory and applications.
  * [6] E. Avila-Vales and Á.G.C. Pérez.  Dynamics of a reaction–diffusion sirs model with general incidence rate in a heterogeneous environment.  Z. Angew. Math. Phys., 73(9), 2022.
  * [7] Arunava Bhadra, Arindam Mukherjee, and Kabita Sarkar.  Impact of population density on covid-19 infected and mortality rate in india.  Modeling Earth Systems and Environment, 7:1–7, 10 2020.
  * [8] Iris Bosa, Adriana Castelli, Michele Castelli, Oriana Ciani, Amelia Compagni, Matteo M. Galizzi, Matteo Garofano, Simone Ghislandi, Margherita Giannoni, Giorgia Marini, and et al.  Response to covid-19: was italy (un)prepared?  Health Economics, Policy and Law, 17(1):1–13, 2022.
  * [9] J. David, A. S. Iyaniwura, M. J. Ward, and F. Brauer.  A novel approach to modelling the spatial spread of airborne diseases: an epidemic model with indirect transmission.  Mathematical Biosciences and Engineering, 17:3294–3328, 2020.
  * [10] Thomas Dierkes, Susanna Röblitz, Moritz Wade, and Peter Deuflhard.  Parameter identification in large kinetic networks with bioparkin, 2013.
  * [11] Patrick E. Farrell, Johan Hake, Simon W. Funke, and Marie E. Rognes.  Automated adjoints of coupled pde-ode systems.  SIAM Journal on Scientific Computing, 41(3), 2019.
  * [12] A. Fernández-Romero, Francisco Guillén-González, and Antonio Suárez.  Theoretical and numerical analysis for a hybrid tumor model with diffusion depending on vasculature.  Journal of Mathematical Analysis and Applications, 503(2), 2021.
  * [13] Félix Foutel-Rodier, François Blanquart, Philibert Courau, Peter Czuppon, Jean-Jil Duchamps, Jasmine Gamblin, Élise Kerdoncuff, Rob Kulathinal, Léo Régnier, Laura Vuduc, et al.  From individual-based epidemic models to mckendrick-von foerster pdes: A guide to modeling and inferring covid-19 dynamics.  Journal of Mathematical Biology, 85(4):43, 2022.
  * [14] Miha Fošnarič, Tina Kamenšek, Jerneja Gros, and Janez Žibert.  Extended compartmental model for modeling covid-19 epidemic in slovenia.  Scientific Reports, 12, 10 2022.
  * [15] Henri P. Gavin.  The levenberg-marquardt algorithm for nonlinear least squares curve-fitting problems.  2019\.
  * [16] Gerrit Großmann, Michael Backenköhler, and Verena Wolf.  Why ode models for covid-19 fail: Heterogeneity shapes epidemic dynamics.  medRxiv, 2021.
  * [17] Nicola Guglielmi, Elisa Iacomini, and Alex Viguerie.  Delay differential equations for the spatially resolved simulation of epidemics with specific application to covid-19.  Mathematical Methods in the Applied Sciences, 45:n/a–n/a, 01 2022.
  * [18] S. Götschel, A. Schiela, and M. Weiser.  Kaskade 7 – a flexible finite element toolbox.  Comp. Math. Appl., 81:444–458, 2020.
  * [19] S. A. Iyaniwura and M. J. Ward.  Synchrony and oscillatory dynamics for a 2-D PDE-ODE model of diffusion-sensing with small signaling compartments.  2020\.
  * [20] Kaitlyn E. Johnson, Grant Howard, William Mo, Michael K. Strasser, Ernesto A. B. F. Lima, Sui Huang, and Amy Brock.  Cancer cell population growth kinetics at low densities deviate from the exponential growth model and suggest an allee effect.  PLOS Biology, 17(8):1–29, 08 2019.
  * [21] Bipin Kumar and Rajesh Kumar Sinha.  Dynamics of an eco-epidemic model with allee effect in prey and disease in predator.  Computational and Mathematical Biophysics, 11(1):20230108, 2023.
  * [22] Y. Kwak, J. Hwang, and C. Yoo.  A new damping strategy of levenberg-marquardt algorithm for multilayer perceptrons.  Neural Network World, 21:327–340, 2011.
  * [23] Thibault Liard, Raphael Stern, and Maria Laura Delle Monache.  A pde-ode model for traffic control with autonomous vehicles.  Networks and Heterogeneous Media, 18(3), 2023.
  * [24] Rachel Mester, Alfonso Landeros, Chris Rackauckas, and Kenneth Lange.  Differential methods for assessing sensitivity in biological models.  PLOS Computational Biology, 18(6):1–30, 06 2022.
  * [25] Sebastian A. Müller, Michael Balmer, William Charlton, Ricardo Ewert, Andreas Neumann, Christian Rakow, Tilmann Schlenther, and Kai Nagel.  Predicting the effects of covid-19 related interventions in urban settings by combining activity-based modelling, agent-based simulation, and mobile phone data.  PLoS One, 16(10), 2021.
  * [26] J.-H. Niemann, S. Uram, S. Wolf, N. Djurdjevac Conrad, and M. Weiser.  Multilevel optimization for policy design with agent-based epidemic models.  J. Comput. Sci., page 102242, 2024.
  * [27] J.R. Rice, E.A. Vavalis, and Daoqi Yang.  Analysis of a nonoverlapping domain decomposition method for elliptic partial differential equations.  Journal of Computational and Applied Mathematics, 87(1):11–19, 1997.
  * [28] G. Schneckenreither, N. Popper, G. Zauner, and F. Breitenecker.  Modelling sir-type epidemics by odes, pdes, difference equations and cellular automata – a comparative study.  Simulation Modelling Practice and Theory, 16(8):1014–1023, 2008.  EUROSIM 2007.
  * [29] Jonathan Richard Shewchuk.  Software triangle.  https://www.cs.cmu.edu/~quake/triangle.html.
  * [30] Sadegh Sabouri Shima Hamidi and Reid Ewing.  Does density aggravate the covid-19 pandemic?  Journal of the American Planning Association, 86(4):495–509, 2020.
  * [31] Software ParaView.  ParaView.  https://www.paraview.org.
  * [32] P.P. Turnes and L.H.A. Monteiro.  An epidemic model to evaluate the homogeneous mixing assumption.  Communications in Nonlinear Science and Numerical Simulation, 19(11):4042–4047, 2014.
  * [33] A. Viguerie, G. Lorenzo, F. Auricchio, D. Baroli, T.J.R. Hughes, A. Patton, A. Reali, T.E. Yankeelov, and A. Veneziani.  Simulating the spread of covid-19 via a spatially-resolved susceptible–exposed–infected–recovered–deceased (seird) model with heterogeneous diffusion.  Applied Mathematics Letters, 111:106617, 2021.
  * [34] Alex Viguerie, Alessandro Veneziani, Guillermo Lorenzo Gómez, Davide Baroli, Nicole Aretz, Alessia Patton, Thomas Yankeelov, Alessandro Reali, Thomas Hughes, and Ferdinando Auricchio.  Diffusion–reaction compartmental models formulated in a continuum mechanics framework: application to covid-19, mathematical analysis, and numerical study.  Computational Mechanics, 66:1–22, 11 2020.
  * [35] David W. S. Wong and Yun Li.  Spreading of covid-19: Density matters.  PLOS ONE, 15(12):1–16, 12 2020.
  * [36] Hanna Wulkow, T. O. F. Conrad, Natasa Djurdjevac Conrad, Sebastian Alexander Mueller, Kai Nagel, and Ch. Schütte.  Prediction of covid-19 spreading and optimal coordination of counter-measures: From microscopic to macroscopic models to pareto fronts.  PLoS ONE, 16(4), 2021.
  * [37] Jinyao Zhao, Xuejuan Zhang, and Jinling Zhao.  A levenberg-marquardt method for tensor approximation.  Symmetry, 15(3), 2023.

Appendices

##  Appendix A Levenberg-Marquardt Algorithm

Algorithm 1 Levenberg-Marquardt algorithm

1:procedure LM($F(p^{(0)}),\text{DF}(p^{(0)}),p^{(0)}$)

2: $l\leftarrow 0$

3: $\lambda_{2}\leftarrow 1$

4: $\lambda_{2}^{\text{positive}}\leftarrow 1$

5: $||F(p^{(-1)})||_{2}\leftarrow-1$

6: repeat

7: choose $\lambda_{1}\geq 0$ so that
$\text{DF}(p^{(l)})^{T}\text{DF}(p^{(l)})+\lambda_{1}I$ is strictly diagonally
dominant

8: $\Delta
p^{(l)}\leftarrow\text{solve}\bigl{(}\text{DF}(p^{(l)})^{T}\text{DF}(p^{(l)})+\lambda_{1}I,-\text{DF}(p^{(l)})^{T}F(p^{(l)})\bigr{)}$

9:

10: decreased_resi $\leftarrow\text{true}$

11: if $||F(p^{(l-1)})||_{2}\geq 0$ then

12: if $||F(p^{(l)})||_{2}<||F(p^{(l-1)})||_{2}$ then

13: $\lambda_{2}\leftarrow 1$ $\triangleright$ reset damping coefficient

14: else if $||F(p^{(l)})||_{2}$ sufficiently small after some runs then

15: $\lambda_{2}\leftarrow 1$ $\triangleright$ reset because $\lambda_{2}$ too
small

16: else

17: decreased_resi $\leftarrow\text{false}$

18: $\lambda_{2}\leftarrow\frac{1}{2}\lambda_{2}$ $\triangleright$ choose
smaller damping coefficient

19: end if

20: end if

21:

22: if not decreased_resi then

23: $l\leftarrow l-1$ $\triangleright$ go back and repeat with new choice of
parameters

24: end if

25:

26: $\lambda_{2}^{\text{positive}}\leftarrow\lambda_{2}$

27: $p^{(l+1)}\leftarrow p^{(l)}+\lambda_{2}^{\text{positive}}\Delta p^{(l)}$

28: while $p^{(l+1)}\leq 0$ do $\triangleright$ in case new parameters are not
positive

29:
$\lambda_{2}^{\text{positive}}\leftarrow\frac{1}{2}\lambda_{2}^{\text{positive}}$

30: $p^{(l+1)}\leftarrow p^{(l)}+\lambda_{2}^{\text{positive}}\Delta p^{(l)}$

31: end while

32:

33: if not decreased_resi then

34: $\lambda_{2}\leftarrow\lambda_{2}^{\text{positive}}$

35: end if

36: $p^{(l)}\leftarrow p^{(l+1)}$ $\triangleright$ continue with new iteration

37: $l\leftarrow l+1$

38: until $||F(p^{(l)})||_{2}$ is sufficiently small or number of maximum
iterations is reached

39:end procedure