1. Introduction

The understanding of disease dynamics for the purpose of prevention and control has become extremely crucial in the recent years. The emergence and reemergence of infectious diseases such as influenza, HIV/AIDS, SARS, and more recently the sudden outburst of the Ebola and Zika virus are events of concern and interest to the general population throughout the world. Moreover, the environmental landscape in which we live is dynamic and often experiences dramatic shifts due to technological innovations that periodically alter the bounds of what we think is possible.

Mathematical models and computer simulations have become irreplaceable experimental tools for building and testing theories, assessing quantitative conjectures, answering specific questions, determining sensitivities to changes in parameter values, and estimating key parameters from data. Understanding the transmission characteristics of infectious diseases in communities, regions, and countries may lead to better approaches to decreasing the transmission of these diseases.

A classical epidemiologal model is the renowned SIR model formulated by Kermack and McKendrick in 1927 [1], [5], [19], [20], [22], [31], which describes the spread of a disease among a single species of N individuals. It is a compartmental model, i.e., the population is split up into three classes of individuals denoted by S, I, R, representing the total number of susceptible, infected and recovered individuals, respectively. Since the effective time period is assumed to be sufficiently short, the model considers neither birth nor death phenomena, nor migration of individuals. It further assumes that susceptible individuals S have never been exposed to the disease, and that they may only be infected by contagious individuals. If \(\beta > 0\) (transmission rate) denotes the average number of adequate contacts of a person per unit time, multiplied by the risk of infection, given contact between an infectious and a susceptible individual, and \(\tau = 1/\gamma\) is the mean waiting time until full recovery, then the SIR model reads

\[\frac{dS}{dt} = -\beta SI, \qquad \frac{dI}{dt} = \beta SI - \gamma I, \qquad \frac{dR}{dt} = \gamma I,\]

supplemented with initial values \(S(0) = S_0\), \(I(0) = I_0\), \(R(0) = R_0\) for some \(S_0\), \(I_0\), \(R_0 \in \mathbb{R}_{\geq 0}\). Clearly the system conserves the number of individuals N = S + I + R for all time \(t \geq 0\). Since its introduction, extensive work has been done on extending the model in various directions.

The understanding of human mobility plays a fundamental role to the research of vector-based and rapid geographical spread of emergent infectious diseases. A popular but rudimentary way to incorporate the

Received August 29^th, 2017, Revised January 8^th, 2018, Accepted for publication January 15^th, 2018. Copyright ©2018 Published by Bio Mathematical Society, e-ISSN: 2549-2896, DOI:10.5614/cbms.2018.1.2.1

spatial movement of hosts into epidemic models is to assume some type of random host movement, leading to reaction-diffusion type equations [32]. This strand of development was built on the pioneering work of Fisher in 1937, who used a logistic-based reaction-diffusion model to investigate the spread of an advantageous gene in a spatially extended population [17]. For considering populations on large geographical scales, scientists have integrated topological features of traffic networks, such as highways, railways and air transportation into models for disease dynamics [25]. Many of these models are stochastic in nature, which results from considering general random walks such as Brownian motion and Lévy flights [7], [25]. In recent years, mathematical methods based on the kinetic theory of active particles has drastically evolved, including the field of mathematical epidemiology [2], [11], [12].

Agent-based models and interacting particle systems have been widely used in understanding how order and stability, or a lack thereof, arises from the interaction of many agents [26]. In addition to the rigorous analysis of models that arise in statistical physics, biology, economics and, now, even in the physics of society, they also provide the means to predict global behavior of a system from the local dynamics between agents [29], [30]. Despite their simplicity, interacting particle systems may be easily extended to include highly complex interactions. Unfortunately, as the number of agents N in the system increases, immense computational cost becomes inevitable.

Multiscale modeling provides a way out. One begins by passing to the so-called mean-field limit \(N \to \infty\), to obtain equations that describe the evolution of the probability density function \(f_t\), over the possible states of the agents [4], [6], [8], [9], [15], [16], [33], [34]. From this probabilistic description, one may further characterize equations corresponding to average/statistical quantities, which have the tendency to be analytically well understood, and numerically tractable.

Fig. 1: Multiscale modeling.

In this paper, we discuss a multiscale approach in deriving macroscopic epidemiology models from agentbased (microscopic) models. Fig. 1 illustrates the general strategy in passing over from the microscopic to macroscopic regime, via the mesoscopic regime, thereby introducing the notion of multiscale [23]. In Section 2 we introduce the agent-based (microscopic) model for epidemiology, that should generalize the classical SIR model in two ways, namely the inclusion of mobility and the continuous transition of an agent's health status. The latter would allow for agents to resist an infection, which is absent from the classical SIR model. A different model for spatially inhomogeneous epidemiology obtained via a multiscale approach has been recently derived in [3]. Section 3 provides an overview and several results pertaining to the limiting equation when the number of agents N tends to infinity. We further provide a way to derive the corresponding macroscopic equation, which is a partial differential equation over the health status. In Section 4 we focus on the well-posedness of the macroscopic equation derived in the previous section and discuss possible stationary distributions. In this section we also point out a way to determine parameters, under which an epidemic may occur for the macroscopic model. Section 5 is devoted to the numerical investigation of the models discussed in the previous sections. Here, we consider various scenarios that justify the adoption of the agent-based epidemiology model and its macroscopic counterpart to model disease dynamics in a spatially inhomogeneous environment. We finally conclude the paper in Section 6 with an outlook to future work and possible extensions of the models introduced within this paper.

2. AN AGENT-BASED EPIDEMIOLOGY MODEL

In this framework, we consider a system of \(N \in \mathbb{N}\) identical agents with position \(X_t^i \in \Omega \subset \mathbb{R}^d\), where \(\Omega\) is a smooth domain, and activity/health status \(U_t^i \in J \subset \mathbb{R}\), \(i=1,\ldots,N\) at time \(t \in [0,\infty)\), satisfying for \(i \in \{1,\ldots,N\}\), the system of stochastic differential equations

\[dX_t^i = \sqrt{2\sigma} \, dW_t^i, \qquad dU_t^i = -\mathcal{H}'(U_t^i) \, dt + F_N(X_t^i, U_t^i, \mathbf{X}_t, \mathbf{U}_t) \, dt, \tag{1}\] with \(\mathbf{X}_t = (X_t^1, \dots, X_t^N)\), \(\mathbf{U}_t = (U_t^1, \dots, U_t^N)\), the standard Wiener process \(W_t^i \in \mathbb{R}^d\) (with normal reflecting boundary conditions whenever \(\Omega\) is bounded), and

\[F_N(X_t^i, U_t^i, \mathbf{X}_t, \mathbf{U}_t) = \frac{1}{N} \sum_{i \neq i} \mathcal{K}(X_t^i, U_t^i, X_t^j, U_t^j),\] where \(\mathcal{H}\colon J\to\mathbb{R}\) is a given potential landscape describing the transition between two activity status, and \(\mathcal{K}\colon S\times S\to\mathbb{R}\) is the force describing inter-agent interactions on the state space \(S:=\Omega\times J\). The initial configuration of the \(N\in\mathbb{N}\) agents in \(S^N\) are assumed to be independent and identically distributed random variables.

Mimicking the SIR model, we set J=[-1,1]. Then the activity \(U^i_t\) is the internal variable describing an agent's health status, where \(U^i_t=-1\) is considered immune due to innate immunity, while \(U^i_t=-1\) is said to be fully recovered, and therefore also immune. This gives us the possibility to work with a continuous health state, which mimics reality, since the infection of an individual intensifies or diminishes more or less continuously. As a first spatial extension to the SIR model and for clarity of presentation, the agents are only modelled to move randomly following the standard Wiener process with diffusion coefficient \(\sqrt{2\sigma}\).

Remark 1. The factor 1/N in \(F_N\) rescales the interaction force K, and is typically called the weak coupling scaling, which will allow for the passage to the mean-field. We will provide examples in Section 5 for this case. Other forms of rescaling may be possible, but will not be discussed here.

The underlying idea to prescribe \(\mathcal{H}\) and \(\mathcal{K}\) stems from reaction rate theory and is adapted to disease dynamics (cf. [18] and references therein). The following are three phenomenological features that are accounted for in our current model:

• 1) A susceptible agent remains susceptible unless exposed to infectious agents. The exposure would need to exceed a certain threshold \(u_* := \arg \max_{u \in J} \mathcal{H}(u)\) for an agent to inherit the disease. This provides the possibility for an agent to resist the disease.
• 2) If a susceptible agent is exposed to a sufficient amount of infectious agents over a duration of time, it will exceed the threshold \(u_*\) and become infected. From that moment on, the recovery phase begins. The agent will be infectious for a period of time and then lose its ability to infect others when its activity exceeds some \(\bar{u} \in (u_*, 1)\). After some time it will arrive at the recovered state.
• 3) Having reached the recovered state, the agent becomes immune to the disease. Therefore u=1 should be an attractor, i.e., local minimum of \(\mathcal{H}\).

For comparison with the classical SIR model, we devide the interval J into compartments, indicating the current active health state. We denote the disjoint partition of J by \(\mathcal{S}=(-1,u_*)\), which represents the susceptible class, \(\mathcal{I}=(u_*,\bar{u})\), the infectious class, and \(\mathcal{R}=(\bar{u},1)\), the fully recovered class. The magnitude of each compartment is then measured simply by counting the number of agents located in the corresponding interval.

Owing to these features mentioned above, we formalize them in the following definition.

Definition 1. A potential landscape \(\mathcal{H}\) is said to be feasible if \(\mathcal{H} \in \operatorname{Lip}_h^1(J)\)

• 1) has two local minima at u = 1 and u = -1, respectively, and
• 2) has one global maximum at \(u_* \in (-1,1)\),

An inter-agent interaction force K is said to be feasible if \(K \in \text{Lip}_b(S \times S)\),

• 1) \(\mathcal{K}(x, u, y, u)\) vanishes for any \(x, y \in \Omega\), \(u \in J\), and
• 2) \(\mathcal{K}(x, u, y, \nu)\) vanishes at \(u \in \mathcal{R} \cup \{-1\}\) for any \(x, y \in \Omega\), \(\nu \in J\).

Remark 2. Notice that an agent at the state u = -1 is considered to be immune to the disease. For this reason, our feasible interaction force \(\mathcal{K}\) does not allow interactions with agents that are at the state u = -1.

The solvability of our microscopic model (1) for any finite number of agents \(N \in \mathbb{N}\) and \(\Omega = \mathbb{R}^d\) is an easy consequence of the strong existence and uniqueness for Itô processes [14]. When \(\Omega\) is bounded, well-posedness may be deduced from eg. [27].

In the following we set \(Z_t^i = (X_t^i, U_t^i)\) for \(1 \le i \le N\) and denote \(\mathcal{P}_p(S)\) to be the set of Borel probability measures with finite p-th moment.

Proposition 1. Let T > 0 be arbitrary, \(\mathcal{H}\) and \(\mathcal{K}\) be feasible, and \(N \in \mathbb{N}\). Furthermore, let the initial values \(\{Z_0^i\}\) be mutually independent and \(f_0\)-distributed random variables with \(f_0 \in \mathcal{P}_2(S)\). Then there exists a unique strong solution \(\mathbf{Z}_t = (\mathbf{X}_t, \mathbf{U}_t)\) of the microscopic model (1) with \(\mathbb{E}\left[\int_0^T |\mathbf{Z}_t|^2 dt\right] < \infty\).

Potential landscapes: An exemplary class of potential landscapes \(\mathcal{H}\) satisfying these features are known as double-well potentials, which includes, for example, potentials of the form

\[\mathcal{H}_{\alpha,\beta}(u) = \alpha(u^2 - 1)^2 + \beta(1 - \sin(\pi u/2)),\] for suitable parameters \(\alpha, \beta \geq 0\). It is easy to see that the set of parameters \(\alpha > 0\), \(\beta < 2^5 \alpha/\pi^2\), provides a feasible set of potential landscapes. One observes that the minima are local attractors. Therefore, an agent at the state u in a neighborhood around \(\{-1,1\}\) will remain there unless perturbed by a sufficient amount of external force. Moreover, an agent is said to have been infected at some point in time if its activity exceeded the threshold \(\mathcal{H}(u_*)\).

There are other possibilities for the potential landscape \(\mathcal{H}\). For instance, one may use a smooth version of a piecewise affine linear function, as seen in Fig. 2. Practically, the potential landscapes provided in Fig. 2 describes a disease that is easily contracted, and requires a relatively long time for complete recovery.

Fig. 2: Examples of potential landscapes \(\mathcal{H}\).

Interaction forces: The inter-agent interactions are described by the force \(\mathcal{K}\), which typically depends on the distance \(r = |X_t^i - X_t^j|\) between agents i and j, and their corresponding activities \(U_t^i\) and \(U_t^j\). A product ansatz of the form

\[\mathcal{K}(x, u, y, \nu) = \Phi(x - y) \psi(u) \chi(\nu), \qquad (x, u), (y, \nu) \in S,\]

may be used, since it easily captures the features mentioned above. Roughly speaking, the function \(\Phi\) should be a non-negative even function, i.e., \(\Phi(r) = \Phi(-r)\), which indicates if two agents are within close proximity for possible interactions to occur. Since the probability of infection is highest when particles are closest, we set \(\Phi(0) = 1\). On the other hand, the function \(\psi\) indicates when an agent is susceptible to infection, whereas \(\chi\) indicates when an agent is infectious. Based on our basic features, an infected agent should not be allowed to infect an agent that is already infected. Furthermore, an infected agent should not change the activity of a recovered agent. Therefore, a feasible \(\mathcal K\) would have that \(\operatorname{supp}(\psi) \cap \operatorname{supp}(\chi) = \emptyset\). A sufficient condition for this would be to require \(\operatorname{supp}(\psi) \subset \mathcal S\) and \(\operatorname{supp}(\chi) \subset \mathcal I\).

To conceive an easy and intuitive interaction force that fits the requirements for the model, we choose \(\Phi(r) = \mathbbm{1}_{B_R(0)}(r)\), \(\psi = \mathbbm{1}_{\mathcal{S}}\) and \(\chi = c_\chi \mathbbm{1}_{\mathcal{I}}\), where \(c_\chi\) represents the strength of infection. Fig. 3 provides an elementary example of an interaction force \(\mathcal{K}\) made up of the functions \(\Phi, \psi\) and \(\chi\). For feasibility reasons, smooth versions of \(\Phi\), \(\psi\) and \(\chi\) shown in the figure are used instead.

3. THE MEAN-FIELD AND MACROSCOPIC EQUATIONS

In this section, we discuss the limiting process that appears when passing to the mean-field limit \(N \to \infty\). The idea in obtaining a limiting equation is to replace the interaction term \(F_N\), which depends on all binary interactions of any pair \((Z_t^i, Z_t^j)\), with an interaction term \(\mathcal{F}\) that describes the interaction of a single agent

Fig. 3: Possible choice of interaction force \(\mathcal{K}(x, u, y, \nu) = \Phi(x - y) \psi(u) \chi(\nu)\).

with an averaged field, the so-called mean-field \(f_t\). Intuitively, if one considers the empirical measure \(\mu_t^N\) of the stochastic processes \(\{Z_t^i\}\) given by

\[\mu^N_t(dz) = \frac{1}{N} \sum\nolimits_{i=1}^N \delta_{Z^i_t}(dz),\] where \(\delta_z\) denotes the Dirac measure at \(z \in S\), then one could formally write

\[F_N(Z_t^i, \mathbf{Z}_t) = \frac{1}{N} \sum_{j \neq i} \mathcal{K}(Z_t^i, Z_t^j) = \int_S \mathcal{K}(Z_t^i, z') \mu_t^N(dz').\]

One then strives to show that the empirical measure \(\mu_t^N\) converges to a deterministic limit measure \(f_t\) in the sense of convergence in law for the underlying random variables. In this case, one may then show that

\[F_N(Z_t^i, \mathbf{Z}_t) \longrightarrow \mathcal{F}[f_t](Z_t^i) = \int_S \mathcal{K}(Z_t^i, z') f_t(dz'),\] in some appropriate notion of convergence, as \(N \to \infty\).

3.1. Nonlinear process and mean-field equations

Indeed, the conjecture is that the process \(Z_t^i\) for some fixed \(1 \le i \le N\) generated by the microscopic system (1) converges to a mean-field process, the so-called McKean nonlinear process, given by the solution of

\[d\bar{X}_t^i = \sqrt{2\sigma} \, dW_t^i, \qquad d\bar{U}_t^i = -\mathcal{H}'(\bar{U}_t^i) \, dt + \mathcal{F}[f_t](\bar{X}_t^i, \bar{U}_t^i) \, dt, \tag{2}\] with the same family of standard Wiener processes \(\{W_t^i\}\) as in (1), and

\[\mathcal{F}[f_t](z) = \int_{S} \mathcal{K}(z, z') f_t(dz'),\] where \(f_t = \text{law}(\bar{Z}^i_t)\) is the law of the random variable \(\bar{Z}^i_t = (\bar{X}^i_t, \bar{U}^i_t)\). Since the law \(f_t\) is required in the definition of the process \(\bar{Z}^i_t\), we have a nonlinear stochastic system at hand. This nonlinear process is supplemented with mutually independent \(f_0\)-distributed initial conditions \(\bar{Z}^i_0\). Note that the solutions \(\{\bar{Z}^i_t\}\) are also independent and identically distributed with the joint law \(f_t^{\otimes N}\).

Applying Itô's formula to the nonlinear process provides an evolution equation for their common law \(f_t\), given by

\[\partial_t f_t - \partial_u (\mathcal{H}' f_t - \mathcal{F}[f_t] f_t) = \sigma \Delta_x f_t, \qquad \lim_{t \to 0} f_t = f_0. \tag{3}\]

This nonlinear and nonlocal kinetic equation is commonly known as the Fokker–Planck equation corresponding to the nonlinear process (2).

As in the microscopic case, we recall an existence and uniqueness result for the nonlinear process (2), as well as the nonlinear kinetic equation (3). Under the feasibility assumptions on \(\mathcal{H}\) and \(\mathcal{K}\), the proof of the following result is rather standard and may be found, for example, in [4], [34].

Proposition 2. Let T > 0 be arbitrary, \(\mathcal{H}\) and \(\mathcal{K}\) be feasible, \(f_0 \in \mathcal{P}_2(S)\). Then the nonlinear process (2) has a pathwise unique solution \(\bar{Z} \in \mathcal{C}([0,T),S)\), with \(f_t = law(\bar{Z}_t) \in \mathcal{P}_2(S)\), \(t \in [0,T)\), satisfying the Fokker-Planck equation (3) in the sense of distributions.

Having unique strong solutions corresponding to (1) and (2), we may provide a quantitative estimate of the difference between the two solutions \(Z^i_t\) and \(\bar{Z}^i_t\) for any \(1 \le i \le N\), and consequently also the difference between their respective laws. For completeness, we provide the proof of the following theorem in Appendix A.

Theorem 1. Let T > 0 be arbitrary, \(\mathcal{H}\) and \(\mathcal{K}\) be feasible, and \(f_0 \in \mathcal{P}_2(S)\). Consider the solutions \(Z_t^i\), \(\bar{Z}_t^i\) to the equations (1), (2) for \(t \in [0,T]\) and each \(1 \le i \le N \in \mathbb{N}\) with mutually independent \(f_0\)-distributed \(Z_0^i\), \(\bar{Z}_0^i\), provided by Propositions 1 and 2 respectively. Then there exists a constant C > 0, independent of \(N \in \mathbb{N}\), such that

\[\sup_{t \in [0,T]} \mathbb{E}[|Z_t^i - \bar{Z}_t^i|^2] \le CN^{-1},\tag{4}\] for any \(1 \le i \le N \in \mathbb{N}\).

The property of the stochastic empirical measure becoming deterministic in the limit is equivalent to the requirement that the law of the N particles become chaotic in the limit [34]. This means that, for a fixed k, the law of the first k agents \(f_t^{(k)}\) satisfies

\[f_t^{(k)} \longrightarrow f_t^{\otimes k} \quad \text{in } \mathcal{P}(S^k),\] as \(N \to \infty\), assuming the k agents to be initially \(f_0^{\otimes k}\)-distributed.

In fact, estimate (4) ensures both theses properties:

1) Propagation of chaos property. Indeed, we deduce from (4) the estimate

\[W_2^2(f_t^{(k)}, f_t^{\otimes k}) \le \mathbb{E}[|(Z_t^1, \dots, Z_t^k) - (\bar{Z}_t^1, \dots, \bar{Z}_t^k)|^2] \le k C N^{-1}\] where \(W_2\) denotes the Wasserstein distance between measures in \(\mathcal{P}_2(S)\) defined by

\[W_2(\mu, \bar{\mu}) = \inf \sqrt{\mathbb{E}[|Z - \bar{Z}|^2]}\]

The infimum is taken over all coupling of random variables \((Z, \bar{Z})\) in \(S \times S\) having distributions \(\mu\) and \(\bar{\mu}\) respectively (cf. [36]).

2) Convergence of the stochastic empirical measure \(\mu_t^N\) towards the deterministic mean-field distribution \(f_t\). Due to (4), we have for any \(\varphi \in \operatorname{Lip}_b(S)\) the estimate

\[\begin{split} \mathbb{E}\left[\left|\frac{1}{N}\sum\nolimits_{i=1}^{N}\varphi(Z_{t}^{i})-\int_{S}\varphi f_{t}(dz)\right|\right] \\ &\leq 2\mathbb{E}\left[\frac{1}{N}\sum\nolimits_{i=1}^{N}|\varphi(Z_{t}^{i})-\varphi(\bar{Z}_{t}^{i})|^{2}+\left|\frac{1}{N}\sum\nolimits_{i=1}^{N}\varphi(\bar{Z}_{t}^{i})-\int_{S}\varphi f_{t}(dz)\right|^{2}\right] \leq CN^{-1}, \end{split}\] for some constant C>0 independent of N and \(t\in[0,T]\). Notice that the second term in the first inequality follows from the law of large numbers. Indeed, this holds since \(\{\bar{Z}_t^i\}\) are mutually independent and identically distributed (see also Appendix A).

Rigorous results of this form were known for the deterministic case since the 70's [6], [13], [33], and then extended to the stochastic case in [4], [34], see also [28].

3.2. Macroscopic equations

At this point, one may derive equations governing macroscopic quantities based on the moments of \(f_t\) by introducing closure relations or further assumptions on \(\mathcal{H}\) and \(\mathcal{K}\). In the following, we assume \(f_t\) to have a sufficiently smooth density with respect to the Lebesgue measure on S, which we denote again by \(f_t\).

The following are several examples that may be of interest:

Model 1. The zeroth order moment of \(f_t\) w.r.t. u, i.e., the first marginal of \(f_t\):

\[\rho_t = \int_I f_t(\cdot, u) \, du,\]

satisfies the simple heat equation

\[\partial_t \rho_t = \sigma \Delta_x \rho_t, \quad \lim_{t \searrow 0} \rho_t = \rho_0,\] which precisely describes the purely diffusive behavior of the nonlinear stochastic process in its first component, namely \(\bar{X}_t\). Indeed, for feasible \(\mathcal{H}\) and \(\mathcal{K}\), we have

\[\int_{J} \partial_{u} (\mathcal{H}' f_{t} - \mathcal{F}[f_{t}] f_{t}) du = (\mathcal{H}' f_{t} - \mathcal{F}[f_{t}] f_{t}) \Big|_{-1}^{1} = 0.\]

If we consider the Wiener process \(W_t\) in a bounded domain \(\Omega \subset \mathbb{R}^d\) with reflective boundary conditions, i.e., we allow the motion of agents only within a bounded region \(\Omega\), we obtain the homogeneous Neumann boundary condition for the heat equation. In this case, the unique equilibrium for this equation is the constant \(\rho_{\text{stat}} \equiv 1/|\Omega|\), i.e., the uniform distribution in the x-variable.

Remark 3. Instead of considering a bounded domain \(\Omega \subset \mathbb{R}^d\) with reflecting boundary conditions for the Wiener process, one could introduce a sufficiently smooth and convex confining potential \(V \colon \mathbb{R}^d \to \mathbb{R}\) with a sufficiently strong growth condition, and additionally \(\int_{\mathbb{R}^d} e^{-V/\sigma} \, dx = 1\). In this case, the mean-field spatial process becomes

\[d\bar{X}_t = -\nabla_x V(\bar{X}_t) dt + \sqrt{2\sigma} dW_t\] and the resulting Fokker-Planck equation reads

\[\partial_t f_t - \partial_u (\mathcal{H}' f_t - \mathcal{F}[f_t] f_t) = \operatorname{div}_x (\sigma \nabla_x f_t + f_t \nabla_x V).\]

As in Model 1, we may take the first marginal of \(f_t\) to obtain

\[\partial_t \rho_t = \operatorname{div}_x (\sigma \nabla_x \rho_t + \rho_t \nabla_x V),\] which is the Fokker–Planck equation corresponding to the Ornstein–Uhlenbeck process. Its unique stationary state is simply given by \(\rho_{\text{stat}} = e^{-V/\sigma}\).

Model 2. Disintegrating the joint probability distribution \(f_t = f_t(x, u)\) into its first marginal \(\rho_t\) and the corresponding conditional distribution \(g_t^x\), i.e., \(f_t(x, u) = g_t^x(u)\rho_t(x)\), and inserting this into the mean-field equation yields

\[\partial_t g_t^x - \partial_u \Big( \mathcal{H}' g_t^x - \mathcal{F}[f_t] g_t^x \Big) = \sigma \Big( \Delta_x g_t^x + \nabla_x \ln \rho_t^2 \cdot \nabla_x g_t^x \Big), \tag{5}\] which is a closed equation for \(g_t^x\), given \(\rho_t\). In fact, if one is given a stationary spatial density \(\rho_{\text{stat}}\) of the population, this can be included directly by simply setting \(\rho_t = \rho_{\text{stat}}\). Notice that this equation is nonlocal in the spatial variable, unless further assumptions are made.

Nevertheless, (5) allows for the computation of \(g_t^x\) for any given spatial density \(\rho_t\), i.e., also those that do not necessarily satisfy the heat equation. Therefore, this macroscopic equation is capable of describing disease dynamics in spatially inhomogeneous populations, where the spatial inhomogeneity is provided by an arbitrary time dependent spatial density \(\rho_t\).

Model 3. A crude approximation to localize the spatial variable in (5) would be to neglect the spatial derivatives on the right-hand side, and to use the product ansatz for the interaction term of the form \(\mathcal{K}(x,u,y,\nu) = \delta_x(y)\Psi(u)\chi(\nu)\), which may be justified in the following sense. Suppose we rescale the spatial variable as \(\tilde{x} \sim \varepsilon x\) and the density as \(\tilde{f}_t \sim f_t(\cdot/\varepsilon)\), i.e., we assume that the spatial domain \(\Omega\) is large in comparison to the range of interaction given by \(\Phi\). Then, we obtain the scaled equation (dropping the tildes)

\[\partial_t g_t^x - \partial_u \Big( \mathcal{H}' g_t^x - \mathcal{F}_{\varepsilon}[f_t] g_t^x \Big) = \varepsilon^2 \sigma \Big( \Delta_x g_t^x + \nabla_x \ln \rho_t^2 \cdot \nabla_x g_t^x \Big),\] where

\[\mathcal{F}_{\varepsilon}[f_t](x,u) = \int_{S} \varepsilon^{-d} \,\Phi(y/\varepsilon) \psi(u) \chi(\nu) f_t(x-y,\nu) \, dy d\nu.\]

Now, if \(\Phi\) has a form of a mollifier, then \(\varepsilon^{-d} \Phi(\cdot/\varepsilon)\) converges towards the Dirac \(\delta_0\) in distribution. By assuming \(f_t\) to be sufficiently smooth, we may formally pass to the limit \(\varepsilon \to 0\) to obtain

\[\partial_t g_t^x - \partial_u \left( \mathcal{H}' g_t^x - \rho_t \, \mathcal{M}[g_t^x] g_t^x \right) = 0, \tag{6}\] with the nonlocal (in the activity variable u) operator

\[\mathcal{M}[g](u) = \psi(u) \int_{I} \chi(\nu)g(\nu) \, d\nu. \tag{7}\]

This localization procedure in the spatial variable provides a pointwise description of the activity, which, from the numerical point of view, is advantageous over the complete mean-field equation, since solving for \(g_t^x\) with respect to the spatial variable \(x \in \Omega\) may be carried out in a pointwise manner, independent of the activity variable \(u \in J\).

Model 4. Assuming further that \(\rho_t \equiv \rho_{\text{stat}} = 1/|\Omega|\), which is the case in spatially homogeneous epidemiology models, a spatially independent model is recovered in the form

\[\partial_t g_t - \partial_u \Big( \mathcal{H}' g_t - \rho_{\text{stat}} \mathcal{M}[g_t] g_t \Big) = 0, \tag{8}\] which describes the probability distribution only in the activity variable \(u \in J\). One can further extract other relevant information, such as the probability of finding agents with a certain activity set \(A \in \mathcal{B}(J)\), simply given by \(\int_A g_t du \in [0,1]\). For example, by choosing \(A = \mathcal{S}\), we recover the probability of finding particles that are susceptible at time \(t \geq 0\).

The last two equations (6), (8) are the simplest of the macroscopic equations. Nevertheless, they sufficiently exhibit important characteristics of a basic epidemiological model modulo the spatial resolution. For this reason, we will study these equations in more detail in the next section.

4. THE NONLOCAL SPATIALLY HOMOGENEOUS MACROSCOPIC EQUATION

In this section, we provide an analytical study of the macroscopic equation

\[\partial_t g_t^x - \partial_u \Big( \mathcal{H}' g_t^x - \bar{\rho}(x) \, \mathcal{M}[g_t^x] g_t^x \Big) = 0, \qquad \lim_{t \searrow 0} g_t^x = g_0^x, \quad x \in \Omega,\] where \(\bar{\rho}\) is a given stationary smooth spatial distribution, and \(\mathcal{M}\) is as given in (7). Since this equation may be solved pointwise in \(x \in \Omega\), we consider the simpler variant

\[\partial_t g_t - \partial_u \Big( \mathcal{H}' g_t - \bar{\rho} \mathcal{M}[g_t] g_t \Big) = 0, \qquad \lim_{t \searrow 0} g_t = g_0,\] (9)

where \(\bar{\rho}\) is simply a constant, \(\bar{\rho} > 0\), and \(g_0\) is an initial distribution of activity.

4.1. Existence and uniqueness

The well-posedness of a nonlocal continuity equation such as (9) may be found, for example, in [10]. Nevertheless, for the convenience of the reader, we provide the principal ideas behind the solvability of the equation.

As in the standard method of characteristics for first order partial differential equations, we may derive the characteristic equation corresponding to the continuity equation (9), which reads

\[\frac{d}{dt}U_t(u) = -\mathcal{H}'(U_t(u)) + \bar{\rho}\mathcal{M}[g_t](U_t(u)), \qquad U_0(u) = u \in J.\] (10)

One recognizes that this equation is again of the form of a nonlinear process since the flow \(U_t\) depends on its law \(g_t\). In fact, if \(U_t\) satisfies the characteristic equation (10), then its law \(g_t\) may be represented by the push-forward of the flow \(U_t\), i.e., \(g_t = U_t \# g_0\), or equivalently

\[\int_{J} \varphi(u) dg_{t} = \int_{J} (\varphi \circ U_{t})(v) dg_{0} \quad \text{for all } \varphi \in \mathcal{C}_{b}(J).\]

Therefore, we define the notion of a Lagrangian solution of (9) with initial data g⁰ ∈ P1(J) as a probability measure g ∈ C([0, T],P1(J)) satisfying the push-forward formula g^t = Ut#g⁰ with the flow U ∈ C([0, T] × J, J) satisfying (10). It is known that Lagrangian solutions and weak measure solutions for (9) coincide (cf. [10]). The main result of this section is the following theorem.

Theorem 2. Let T > 0 be arbitrary, g⁰ ∈ P1(J), and H be feasible and ψ, χ ∈ Lip^b (S). Then, there exists a unique Lagrangian solution g ∈ C([0, T],P1(J)) to the equation (9).

Its proof relies on the use of the well-known Banach fixed point theorem [37] for complete metric spaces. For this reason, we consider the space C([0, T],P1(J)), endowed with the distance

\[d(\mu, \nu) = \sup_{t \in [0,T]} W_1(\mu_t, \nu_t), \qquad \mu, \nu \in \mathcal{C}([0,T], \mathcal{P}_1(J)),\] where W¹ denotes the 1-Wasserstein distance, given by

\[W_1(\mu,\nu) = \inf_{\pi \in \Pi(\mu,\nu)} \iint_{J \times J} |u - v| d\pi(u,v), \qquad \mu,\nu \in \mathcal{P}_1(J).\]

Here Π(µ, ν) denotes the set of all measures π with marginals π(·, J) = µ and π(J, ·) = ν. It is known that the 1-Wasserstein distance metrizes the narrow convergence in P1(J), i.e., the convergence tested against continuous and bounded functions, which makes (P1(J), W1) a separable complete metric space, since J is complete [36]. Consequently, the function space C([0, T],P1(J)) endowed with the distance d above is also a separable complete metric space.

Now consider, for any given gˆ ∈ C([0, T],P1(J)), the auxiliary problem

\[\frac{d}{dt}U_t(u) = -\mathcal{H}'(U_t(u)) + \bar{\rho}\mathcal{M}[\hat{g}_t](U_t(u)), \qquad U_0(u) = u \in J.\] (11)

It is easy to see that the right-hand side of the equation is continuous in the temporal variable, and globally Lipschitz-continuous in the activity variable for any feasible functions H, ψ and χ. Therefore, the Picard– Lindelof theorem ¨ , or similarly, the Cauchy–Lipschitz theorem, provides a unique global solution U·(u) ∈ C([0, T], J) for any u ∈ J, and thereby a flow U ∈ C([0, T] × J, J). We then construct a new probability measure g ∈ C([0, T],P1(J)) by means of push-forward, i.e., g^t = Ut#g0, where g⁰ ≡ gˆ⁰ ∈ P1(J).

Consequently, this induces a mapping T : C([0, T],P1(J)) → C([0, T],P1(J)), gˆ 7→ g, which we show to admit a fixed point satisfying the nonlocal continuity equation (9). Before proceeding with the proof of Theorem 2, we provide a stability estimate that will assist in showing the required contracting property of the mapping T .

Lemma 1. Let g, ˆ hˆ ∈ C([0, T],P1(J)) be given and g, h be Lagrangian solutions to

\[\partial_t g_t - \partial_u \Big( \mathcal{H}' g_t - \bar{\rho} \, \mathcal{M}[\hat{g}_t] g_t \Big) = 0, \qquad \partial_t h_t - \partial_u \Big( \mathcal{H}' h_t - \bar{\rho} \, \mathcal{M}[\hat{h}_t] h_t \Big) = 0,\] with initial conditions g⁰ = ˆg⁰ and h⁰ = hˆ ⁰ in P1(J), respectively. Then the estimate

\[W_1(g_t, h_t) \le \left(W_1(g_0, h_0) + c_2 \int_0^t W_1(\hat{g}_s, \hat{h}_s) \, ds\right) e^{c_1 t}\] for all \(t \ge 0\), holds true with positive constants c1, c2, depending only on ρ¯, H, ψ and χ.

Proof: We first note that g^t = Ut#g⁰ and h^t = Vt#h0, where U, V ∈ C([0, T] × J, J) satisfy

\[U_t(u) = u - \int_0^t \mathcal{H}'(U_s(u)) - \bar{\rho} \mathcal{M}[\hat{g}_s](U_s(u)) ds,\]
\[V_t(u) = v - \int_0^t \mathcal{H}'(V_s(v)) - \bar{\rho} \mathcal{M}[\hat{h}_s](V_s(v)) ds,\]

respectively. Now let \(\pi_0 \in \Pi(g_0, h_0)\) be an optimal coupling of \(g_0\) and \(h_0\), and \(\pi_t = (U_t, V_t) \# \pi_0\). Then \(\pi_t \in \Pi(U_t \# g_0, V_t \# h_0) = \Pi(g_t, h_t)\), which is not necessarily optimal. For \(\pi_t\), we have

\[\begin{split} W_{1}(g_{t},h_{t}) &\leq \iint_{J\times J} |u-v| \, d\pi_{t}(u,v) = \iint_{J\times J} |U_{t}(u)-V_{t}(v)| \, d\pi_{0}(\hat{u},\hat{v}) \\ &\leq \iint_{J\times J} |u-v| \, d\pi_{0}(u,v) + \int_{0}^{t} \iint_{J\times J} |\mathcal{H}'(U_{s}(u))-\mathcal{H}'(V_{s}(v))| \, d\pi_{0}(u,v) \, ds \\ &+ \bar{\rho} \int_{0}^{t} \iint_{J\times J} |\mathcal{M}[\hat{g}_{s}](U_{s}(u)) - \mathcal{M}[\hat{h}_{s}](V_{s}(v))| \, d\pi_{0}(u,v) \, ds \\ &= W_{1}(h_{0},g_{0}) + I_{1} + I_{2}. \end{split}\]

To estimate \(I_1\), we simply use the Lipschitz-continuity of \(\mathcal{H}'\) to obtain

\[I_1 \leq L_{\mathcal{H}'} \int_0^t \iint_{J \times J} |U_s(u) - V_s(v)| d\pi_0(u, v) ds = L_{\mathcal{H}'} \int_0^t \iint_{J \times J} |u - v| d\pi_s(u, v) ds.\]

Similarly, we use the Lipschitz-continuity of \(\psi\) to obtain

\[I_{2} = \bar{\rho} \int_{0}^{t} \iint_{J \times J} \left| \psi(U_{s}(u)) \left( \int_{J} \chi(\hat{u}) \, d\hat{g}_{s}(\hat{u}) \right) - \psi(V_{s}(v)) \left( \int_{J} \chi(\hat{v}) \, d\hat{h}_{s}(\hat{v}) \right) \right| \, d\pi_{0}(u, v) \, ds\] \[\leq \bar{\rho} L_{\psi} \|\chi\|_{\infty} \int_{0}^{t} \iint_{J \times J} |U_{s}(u) - V_{s}(v)| \, d\pi_{0}(u, v) \, ds\] \[+ \bar{\rho} \|\psi\|_{\infty} \int_{0}^{t} \left| \int_{J} \chi(\hat{u}) \, d\hat{g}_{s}(\hat{u}) - \int_{J} \chi(\hat{v}) \, d\hat{h}_{s}(\hat{v}) \right| \, ds,\] where we used the fact that \(\psi\) and \(\chi\) are bounded, and \(\hat{g}_t\) and \(\pi\) are a probability measures over J and \(J \times J\), respectively. Concerning the last term, we estimate further to obtain

\[\left| \int_{J} \chi(\hat{u}) \, d\hat{g}_{s}(\hat{u}) - \int_{J} \chi(\hat{v}) \, d\hat{h}_{s}(\hat{v}) \right| \leq \iint_{J \times J} |\chi(\hat{u}) - \chi(\hat{v})| \, d(\hat{g}_{s} \otimes \hat{h}_{s})(\hat{u}, \hat{v})\] \[\leq L_{\chi} \iint_{J \times J} |\hat{u} - \hat{v}| \, d(\hat{g}_{s} \otimes \hat{h}_{s})(\hat{u}, \hat{v}).\]

Putting all the terms together yields

\[W_{1}(g_{t}, h_{t}) \leq W_{1}(h_{0}, g_{0}) + \left(L_{\mathcal{H}'} + \bar{\rho}L_{\psi} \|\chi\|_{\infty}\right) \int_{0}^{t} \iint_{J \times J} |u - v| d\pi_{s}(u, v) ds + \bar{\rho}L_{\chi} \|\psi\|_{\infty} \int_{0}^{t} \iint_{J \times J} |u - v| d(\hat{g}_{s} \otimes \hat{h}_{s})(u, v) ds.\]

Optimizing the right-hand side over all possible couplings in \(\Pi(g_s, h_s)\) and \(\Pi(\hat{g}_s, \hat{g}_s)\) gives

\[W_1(g_t, h_t) \le W_1(g_0, h_0) + c_1 \int_0^t W_1(g_s, h_s) \, ds + c_2 \int_0^t W_1(\hat{g}_s, \hat{h}_s) \, ds,\] with \(c_1 = L_{\mathcal{H}'} + \bar{\rho} L_{\psi} \|\chi\|_{\infty}\) and \(c_2 = \bar{\rho} L_{\chi} \|\psi\|_{\infty}\). From Gronwall's inequality, we finally obtain

\[W_1(g_t, h_t) \le \left(W_1(g_0, h_0) + c_2 \int_0^t W_1(\hat{g}_s, \hat{h}_s) \, ds\right) e^{c_1 t},\] which completes the proof.

We now have all the ingredients necessary to complete the proof of Theorem 2.

Proof of Theorem 2: We consider the mapping \(\mathcal{T}: \mathcal{C}([0,T],\mathcal{P}_1(J)) \to \mathcal{C}([0,T],\mathcal{P}_1(J))\) as discussed above. However, we consider a weighted metric of the form

\[d_{\lambda}(g,h) = \sup_{t \in [0,T]} e^{-\lambda t} W_1(g_t, h_t),\] which is clearly equivalent to the usual metric d, for any \(\lambda > 0\). Therefore, the space \(\mathcal{C}([0,T],\mathcal{P}_1(J))\) endowed with the metric \(d_{\lambda}\) is again a separable complete metric space.

Now let \(g = \mathcal{T}(\hat{g})\) and \(h = \mathcal{T}(\hat{h})\), with \(g_0 = h_0\). Then Lemma 1 provides the estimate

\[W_1(g_t, h_t) = c_2 e^{c_1 t} \int_0^t W_1(\hat{g}_s, \hat{h}_s) \, ds \le (c_2 e^{c_1 t} / \lambda) (e^{\lambda t} - 1) d_{\lambda}(\hat{g}, \hat{h}).\]

Multiplying both sides by \(\exp(-\lambda t)\) and taking the supremum over time \(t \in [0,T]\) yields

\[d_{\lambda}(\mathcal{T}(\hat{q}), \mathcal{T}(\hat{h})) < (c_2 e^{c_1 T} / \lambda) d_{\lambda}(\hat{q}, \hat{h}).\]

Therefore, choosing \(\lambda > c_2 \exp(c_1 T)\) makes \(\mathcal{T}\) a contraction mapping with respect to the metric \(d_{\lambda}\). Finally, we invoke the Banach fixed point theorem to obtain a unique fixed point in the space \(\mathcal{C}([0,T],\mathcal{P}_1(J))\), which satisfies the nonlocal macroscopic equation (9).

In fact, one can further show that if the initial measure \(g_0 \in \mathcal{P}_1^{\text{ent}}(J)\), where \(\mathcal{P}_1^{\text{ent}}(J)\) denotes the space of probability measures with finite first moment that are, additionally, absolutely continuous with respect to the Lebesgue measure and have finite entropy

\[0 \leq \operatorname{Ent}(g) = \int_{I} \left( g \log(g) - g + 1 \right) du < \infty, \qquad g \in \mathcal{P}_{1}^{\operatorname{ent}}(J),\] then \(g_t \in \mathcal{P}_1^{\text{ent}}(J)\) for all times \(t \geq 0\). Indeed, assuming \(\mathcal{H}\), \(\psi\) and \(\chi\) to be feasible, then

\[\begin{split} \frac{d}{dt} \mathrm{Ent}(g_t) &= -\int_J \partial_u g_t \Big( \mathcal{H}'(u) - \bar{\rho} \mathcal{M}[g_t] \Big) du = \int_J g_t \Big( \mathcal{H}''(u) - \bar{\rho} \partial_u \mathcal{M}[g_t] \Big) du \\ &\leq \Big( 1 + \bar{\rho} \|\chi\|_{\infty} \Big) \mathrm{Ent}(g_t) + \int_J \Big( e^{\mathcal{H}''(u)} - 1 \Big) du + \bar{\rho} \|\chi\|_{\infty} \int_J \Big( e^{\psi'(u)} - 1 \Big) du \\ &= c_1 \mathrm{Ent}(g_t) + c_2, \end{split}\] where we integrated by parts in the first two equalities, and applied Young's inequality of the form \(ab \le e^a + b \ln(b) - b\) for \(a, b \in \mathbb{R}\), \(b \ge 0\) in the inequality. Equivalently, we have in integral form

\[\operatorname{Ent}(g_t) \le (\operatorname{Ent}(g_0) + c_2 t) + c_1 \int_0^t \operatorname{Ent}(g_s) \, ds.\]

A simple application of the Gronwall inequality leads to the estimate

\[\operatorname{Ent}(g_t) \le (\operatorname{Ent}(g_0) + c_2 t) e^{c_1 t},\] which shows that \(g_t \in \mathcal{P}_1^{\text{ent}}(J)\) for all times \(t \geq 0\) as asserted. Summarizing, we have

Proposition 3. Let \(\mathcal{H}\) be feasible, \(\psi\), \(\chi \in Lip_b(S)\) and \(g_0 \in \mathcal{P}_1^{ent}(J)\). Then \(g_t \in \mathcal{P}_1^{ent}(J)\) for all times \(t \geq 0\).

4.2. Stationary measures and transitions

Here, we would like to explore the possible stationary states of the nonlocal macroscopic equation (9) and provide an expression similar to the classical SIR model in order to determine the occurrence of an epidemic, or otherwise.

As noted in Remark 2, any agent that begins with the state u=-1 remains there for all times. Therefore, we expect \(\delta_{-1}\) to be a natural stationary measure for (9). In fact, it is not difficult to see that, if \(\mathcal{H}\) is feasible and \(\operatorname{supp}(\psi) \cap \operatorname{supp}(\chi) = \emptyset\), then \(\delta_{u_*}\) and \(\delta_1\) are also stationary measures. Indeed, since every stationary measure should satisfy

\[\int_J \left(\mathcal{H}'(u) - \bar{\rho} \mathcal{M}[g_\infty]\right) \partial_u \varphi(u) \, dg_\infty = 0, \qquad \text{for all } \varphi \in \mathcal{C}_b(J),\]

we simply substitute \(g_{\infty} = \delta_{\sigma}\), \(\sigma \in \{-1, u_*, 1\}\) into the equation and use the fact that \(\mathcal{H}'(\sigma) = 0\), supp\((\psi) \cap \text{supp}(\chi) = \emptyset\), to verify its stationarity.

Definition 2 (Asymptotic stable). A stationary state \(g_{\infty} \in \mathcal{P}_1(J)\) is said to be asymptotically stable if there exists \(\varepsilon > 0\) such that the solution \(g^{\eta}\) of (9) with initial condition \(\eta \in \mathcal{P}(J)\) with \(W_{\infty}(g_{\infty}, \eta) < \varepsilon\) converges narrowly to \(g_{\infty}\) as \(t \to \infty\), where \(W_{\infty}\) is the infinite Wasserstein distance defined by

\[W_{\infty}(\mu,\nu) = \inf\nolimits_{\pi \in \Pi(\mu,\nu)} \Big\{ \pi \text{-esssup}_{(u,v) \in J \times J} |u-v| \Big\}, \qquad \mu,\nu \in \mathcal{P}(J).\]

The following result classifies all asymptotically stable states whenever \(\mathcal{H}\) attains strict local minima at u=-1 and u=1.

Theorem 3. Let \(\mathcal{H}\) be feasible, where \(\mathcal{H}\) has strict local minima at \(u \in \{-1,1\}\), a global maxima at \(u_* \in (-1,1)\), and \(\psi\), \(\chi \in Lip_b(S)\) with \(supp(\psi) \cap supp(\chi) = \emptyset\). Then only measure of the form

\[g_{\infty} = (1 - \alpha)\delta_{-1} + \alpha\delta_1, \qquad \alpha \in [0, 1],\] are asymptotically stable states of the nonlocal macroscopic equation (9).

Proof: Let \(g_{\infty} = (1 - \alpha)\delta_{-1} + \alpha\delta_1\) for some \(\alpha \in [0, 1]\). We show that \(g_{\infty}\) is asymptotically stable. In order to determine \(\varepsilon > 0\) appropriately, we first establish non-empty neighborhoods \(B_{-1} \subset \mathcal{S}\) and \(B_1 \subset \mathcal{R}\) around u = -1 and u = 1, respectively, where \(\mathcal{H}\) is strictly convex. Such neighborhoods exists since \(\mathcal{H}\) has strict local minima at \(u \in \{-1, 1\}\). Moreover, these sets are disjoint since \(\mathcal{S}\) and \(\mathcal{R}\) are disjoint. We then choose

\[0 < \varepsilon < \min\{\operatorname{diam}(B_{-1}), \operatorname{diam}(B_1)\},\\] where the diameter of a set A is given by \(diam(A) = \sup\{|x - y|, x, y \in A\}\).

Suppose for the moment that \(\eta \in \mathcal{P}(J)\) is atomic, i.e., \(\eta = \sum_{u \in I_{\eta}} c_u \delta_u\) with \(\sum_{u \in I_{\eta}} c_u = 1\) for some discrete set \(I_{\eta} \subset J\). Now consider two sets \(N_{-1}\) and \(N_1\). We allocate each \(u \in I_{\eta}\) to one of these sets as follows,

\[u \in N_i \Leftrightarrow |u - j| = \min\{|u + 1|, |u - 1|\}, j \in \{-1, 1\}.\]

In this case, it is not difficult to see that any \(\pi \in \Pi(\eta, g_{\infty})\) satisfies

\[\pi\text{-esssup}_{(u,v) \in J \times J} |u-v| = \pi\text{-esssup}_{(u,v) \in I_n \times \{-1,1\}} |u-v| = \max_{j \in \{-1,1\}} \sup_{u \in N_j} |u-j|.\]

Therefore, the requirement \(W_{\infty}(\eta, g_{\infty}) < \varepsilon\) dictates that

\[\max_{j \in \{-1,1\}} \sup_{u \in N_i} |u - j| < \infty.\]

i.e., \(N_j \subset B_j\), \(j \in \{-1,1\}\), and hence \(supp(\eta) \subset B_{-1} \cup B_1\).

Now consider the Lagrangian solution \(g^{\eta}\) corresponding to (9) with initial condition \(\eta\), or equivalently,

\[\frac{d}{dt}U_t^{\eta}(u) = -\mathcal{H}'(U_t^{\eta}(u)) + \bar{\rho}\mathcal{M}[g_t^{\eta}](U_t^{\eta}(u))\] with \(u \in \text{supp}(\eta)\). Since S and R are disjoint sets, we may consider first \(u \in \text{supp}(\eta) \cap S \subset B_{-1}\). In this case, we have that

\[\frac{d}{dt}U_t^{\eta}(u)\Big|_{t=0} = -\mathcal{H}'(u) + \bar{\rho}\mathcal{M}[\eta](u) = -\mathcal{H}'(u) < 0,\] which says that \(U_t^{\eta}(u)\) remains in \(B_{-1}\) for sufficiently small t>0, due to continuity. Analogously, we can show that \(U_t^{\eta}(u) \in B_1\) for any \(u \in \operatorname{supp}(\eta) \cap \mathcal{R} \subset B_1\) when t>0 is sufficiently small. Therefore, the support of \(g_t^{\eta}\) is contained within \(B_{-1} \cup B_1\) for t>0 sufficiently small. By iterating this argument along the flow \(U_t^{\eta}\), we have that \(\operatorname{supp}(g_t^{\eta}) \subset B_{-1} \cup B_1\) for all times \(t \geq 0\). Equivalently, we have that \(U_t^{\eta}(u) \in B_1 \cup B_{-1}\) for all \(t \geq 0\), for any \(u \in \operatorname{supp}(\eta)\).

The previous discussion implies that \(\mathcal{M}[g_t^{\eta}](v) = 0\) for any \(v \in J\), and hence

\[\frac{d}{dt}U_t^{\eta}(u) = -\mathcal{H}'(U_t^{\eta}(u)), \qquad u \in \operatorname{supp}(\eta).\]

Taking the time derivative of \(\mathcal{H}\) along the flow \(U_t(u)\) yields

\[\frac{d}{dt}\mathcal{H}(U_t^{\eta}(u)) = \mathcal{H}'(U_t^{\eta}(u))\frac{d}{dt}U_t^{\eta}(u) = -|\mathcal{H}'(U_t^{\eta}(u))|^2 < 0 \qquad \text{for all } t \ge 0,\] which says that the flow \(U_t^{\eta}(u)\) minimizes \(\mathcal{H}\) with time. Since \(\mathcal{H}\) is strictly convex in \(B_{-1} \cup B_1\), we have that \(\mathcal{H}'(v) \neq 0\) for any \(v \in B_{-1} \cup B_1\), \(v \notin \{-1,1\}\). Thus,

\[U_t^{\eta}(u) \longrightarrow \begin{cases} -1, & \text{for } u \in \text{supp}(\eta) \cap \mathcal{S} \\ 1, & \text{for } u \in \text{supp}(\eta) \cap \mathcal{R} \end{cases} \quad \text{as } t \to \infty.\]

Since \(\operatorname{supp}(\eta) \cap \mathcal{S}\) and \(\operatorname{supp}(\eta) \cap \mathcal{R}\) are disjoint, the mass of \(g_t^{\eta}\) within \(B_{-1}\) is \((1 - \alpha)\), and \(B_1\) is \(\alpha\), due to conservation of mass. Consequently,

\[\int_{J} \phi(u) \, dg_{t}^{\eta} = \int_{J} \phi(u) \, dU_{t}^{\eta} \# \eta = \int_{J} \phi \circ U_{t}^{\eta}(u) \, d\eta\]\[= \int_{B_{-1}} \phi \circ U_{t}^{\eta}(u) \, d\eta + \int_{B_{1}} \phi \circ U_{t}^{\eta}(u) \, d\eta,\] for every \(\phi \in \mathcal{C}_b(J)\). Passing to the limit \(t \to \infty\), we obtain

\[\lim_{t \to \infty} \int_J \phi(u) \, dg_t^{\eta} = \int_{B_{-1}} \phi(-1) \, d\eta + \int_{B_1} \phi(1) \, d\eta = (1 - \alpha)\phi(-1) + \alpha\phi(1) = \int_J \phi \, dg_{\infty},\]

i.e., \(g_t^{\eta}\) converges narrowly to \(g_{\infty}\) as \(t \to \infty\).

Now let \(g_{\infty}\) be an asymptotically stable state. Then \(g_{\infty}\) takes the form

\[g_{\infty} = c_{-1}\delta_{-1} + c_{u_*}\delta_{u_*} + c_1\delta_1\] with \(\sum_{i \in \{-1, u_*, 1\}} c_i = 1\), since these form all possible stationary states. We show that \(c_{u_*} = 0\). To do so, we suppose otherwise, i.e., \(c_{u_*} > 0\), and show that

\[\forall \varepsilon > 0, \exists \eta \text{ with } W_{\infty}(\eta, g_{\infty}) < \varepsilon : \quad g_t^{\eta} \not\to g_{\infty} \quad \text{narrowly as } t \to \infty,\]

i.e., \(g_{\infty}\) is not asymptotically stable. More precisely, we show that

\[\exists q>0, \exists \phi\in\mathcal{C}_b(J): \quad \left|\int_I \phi\,dg_t^\eta - \int_I \phi\,dg_\infty\right| \geq q \ \text{ for all } t\geq 0.\]

We begin by taking an arbitrary \(\varepsilon > 0\). Choosing \(\eta = c_{-1}\delta_{-1} + c_{u_*}\delta_v + c_1\delta_1\), we find that

\[W_{\infty}(\eta, g_{\infty}) \le |v - u_*|.\]

Since \(u^*\) is a global maximum, there is a neighborhood \(B_*\) around \(u_*\) such that \(\mathcal{H}|_{B_*}\) is concave. Therefore, taking \(v \in B_* \cap \mathcal{I}\) such that \(0 < q := |v - u_*| < \varepsilon\) provides a candidate. Now consider the Lagrangian solution \(g_t^\eta = U_t^\eta \# \eta\) where

\[\frac{d}{dt}U_t^{\eta}(u) = -\mathcal{H}'(U_t^{\eta}(u)) + \bar{\rho}\mathcal{M}[g_t^{\eta}](U_t^{\eta}(u)).\]

Following a similar argument as above, we determine that

\[\frac{d}{dt}U_t^{\eta}(v)\Big|_{t=0} = -\mathcal{H}'(v) + \bar{\rho}\mathcal{M}[\eta](v) = -\mathcal{H}'(v) \ge 0,\]

thus implying the existence of some \(t_*>0\) such that \(U^\eta_t(v)\geq v\) for \(t\in[0,t_*)\). Iterating this procedure, we find that \(U^\eta_t(v)\geq v\) for all \(t\geq 0\). Finally, choosing \(\phi(u)=u\), we obtain

\[\left| \int_{I} \phi \, dg_{t}^{\eta} - \int_{I} \phi \, dg_{\infty} \right| = c_{u_{*}} |U_{t}^{\eta}(v) - u_{*}| \ge c_{u_{*}} |v - u_{*}| = q > 0,\] for all \(t \ge 0\), which shows that \(g_t^{\eta} \not\to g_{\infty}\) narrowly as \(t \to \infty\), thus concluding the proof.

Remark 4. One easily verifies that stable stationary states for the mean-field equation (3) may be identified with the measure \(f_{\infty}(dx,du) = \rho_{\text{stat}}(dx) \otimes g_{\infty}(du)\).

We now proceed to derive an equivalent expression for the basic reproduction number \(\mathfrak{R}_0\) present in the classical SIR model, which determines if a disease leads to an epidemic or otherwise. For the classical

SIR model, the basic reproduction number \(\Re_0\) is given by the formula \(\Re_0 = \beta S_0/\gamma\), where \(\beta > 0\) is the transmission rate, \(\gamma\) the recovery rate, and \(S_0\) the initial susceptible population.

To provide a correspondence between the nonlocal macroscopic model (9) and the SIR model, we make simplifying assumptions on \(\mathcal{H}\), \(\psi\) and \(\chi\). More precisely, we assume that

\[\psi = \mathbb{1}_{\mathcal{S}}^{\varepsilon}, \qquad \chi = c_{\chi} \mathbb{1}_{\mathcal{I}}^{\varepsilon}, \qquad \mathcal{H}' = \lambda \mathbb{1}_{\mathcal{S}}^{\varepsilon} - \gamma \mathbb{1}_{\mathcal{I}}^{\varepsilon}, \tag{12}\] where \(c_{\chi}\), \(\lambda\), \(\gamma\) are positive constants, and \(\mathbb{1}_{A}^{\varepsilon}\) are mollified versions of the indicator function over a Borel set \(A \subset J\) with supp\((\mathbb{1}_{A}^{\varepsilon}) \subset A\). We further assume that \(\bar{u} = 1\), i.e., \(\mathcal{R} = \{1\}\).

Definition 3. We define the effective transition from the class of susceptible agents S to the class of infectious agents I as

\[\mathcal{E}_t = \int_{\mathcal{S}} u \, dg_t - \int_{\mathcal{I}} u \, dg_t,\] for all times \(t \ge 0\). Roughly speaking, \(\mathcal{E}_t\) gives an indication of the probability of agents that lie within a infinitesimal neighborhood of \(u_*\), i.e., around the point of transition.

Similar to the classical case, the disease is said to be epidemic if

\[\left. \frac{d}{dt} \mathcal{E}_t \right|_{t=0} > 0,\] which suggests the presence of agents transitioning from class S to class I.

By taking the temporal derivative of \(\mathcal{E}_t\), we obtain

\[\begin{split} \frac{d}{dt}\mathcal{E}_{t} &= -\int_{\mathcal{S}} \left( \mathcal{H}' - \bar{\rho} \,\mathcal{M}[g_{t}] \right) dg_{t} + u \Big( \mathcal{H}'g_{t} - \bar{\rho} \,\mathcal{M}[g_{t}]g_{t} \Big) \Big|_{-1}^{u_{*}} \\ &+ \int_{\mathcal{I}} \left( \mathcal{H}' - \bar{\rho} \,\mathcal{M}[g_{t}] \right) dg_{t} - u \Big( \mathcal{H}'g_{t} - \bar{\rho} \,\mathcal{M}[g_{t}]g_{t} \Big) \Big|_{u_{*}}^{1} \\ &= -\lambda S_{t} + \bar{\rho} \, c_{\chi} S_{t} I_{t} - \gamma I_{t}, \end{split}\] where we denote \(S_t = \int \mathbb{1}_{S}^{\varepsilon} dg_t\) and \(I_t = \int \mathbb{1}_{T}^{\varepsilon} dg_t\). Consequently, we have

\[\frac{d}{dt}\mathcal{E}_t\Big|_{t=0} = \gamma I_0(\mathfrak{R}_0 - 1),\] with the basic reproduction number \(\mathfrak{R}_0 = \bar{\rho} \, c_{\chi} S_0 / \gamma - \lambda S_0 / (\gamma I_0)\), which indicates that an epidemic only occurs when \(\mathfrak{R}_0 > 1\). Notice that if \(\lambda = 0\), we recover the classical basic reproduction number.

Summarizing the discussion above yields the following statement.

Theorem 4. Let \(\mathcal{H}\), \(\psi\) and \(\chi\) be feasible and satisfy additionally (12). Then

\[\frac{d}{dt}\mathcal{E}_t\Big|_{t=0} = \gamma I_0(\mathfrak{R}_0 - 1),\] with the basic reproduction number \(\mathfrak{R}_0 = \bar{\rho} c_{\chi} S_0 / \gamma - \lambda S_0 / (\gamma I_0)\). In particular, an epidemic occurs when \(\mathfrak{R}_0 > 1\).

5. Numerical Investigations

We recall the two epidemiological models that will be under investigation within this section, namely the microscopic model

\[dX_t^i = \sqrt{2\sigma} \, dW_t^i, \qquad dU_t^i = -\mathcal{H}'(U_t^i) \, dt + F_N(X_t^i, U_t^i, \mathbf{X}_t, \mathbf{U}_t) \, dt,\] and the macroscopic model

\[\partial_t g_t - \partial_u \Big( \mathcal{H}' g_t - \bar{\rho} \, \mathcal{M}[g_t] g_t \Big) = 0.\]

In all our numerical simulations, we consider the spatial dynamics to be within the bounded domain \(\Omega = [0,1]^2\), with reflecting boundary conditions for the Wiener processes \(W_t^i\). Furthermore, we only consider interactions of product form, i.e.,

\[\mathcal{K}(x, u, y, \nu) = \Phi(x - y)\psi(u)\chi(\nu), \qquad (x, u), (y, \nu) \in S.\]

The standard Euler-Maruyama scheme was employed to solve the microscopic equations numerically [24]. Appropriate step sizes were chosen to ensure stability of the explicit scheme. In the absence of noise, i.e., \(\sigma=0\), we simply use the standard \(explicit\ Euler\) scheme. Since stochastic processes admit different solution paths for different realizations, we consider multiple realizations (often M=100 realizations) to obtain statistical information such as the mean and variance.

To compute the mean \(m_t^{\varphi}\) of a given observable bounded \(\varphi\), we use the well-known estimator

\[m_t^{\varphi} = \frac{1}{M} \sum_{j=1}^{M} \varphi(\mathbf{Z}_t^{(j)}),\] which ensures convergence towards the mean via the law of large numbers. Here, the superscript index j represents the jth realization of the microscopic simulation. Typical observables we often use are the number of agents within the health classes \(\mathcal{S}\), \(\mathcal{I}\) and \(\mathcal{R}\):

\[S_t = \sum_{i=1}^N \mathbb{1}_{\mathcal{S}}(U_t^i), \qquad I_t = \sum_{i=1}^N \mathbb{1}_{\mathcal{I}}(U_t^i), \qquad R_t = \sum_{i=1}^N \mathbb{1}_{\mathcal{R}}(U_t^i).\]

As an estimator for the variance, we choose the unbiased sample variance

\[\operatorname{var}_t^{\varphi} = \frac{1}{M-1} \sum\nolimits_{i=1}^{M} \left( \varphi(\mathbf{Z}_t^{(j)}) - m_t^{\varphi} \right)^2, \qquad s_t^{\varphi} = \sqrt{\operatorname{var}_t^{\varphi}}.\]

In all the plots below, we use the color blue to identify \(S_t\), red for \(I_t\) and green for \(R_t\). The standard deviation for each observable \(s_t^{\varphi}\) will be shown as shaded regions around its sample mean \(m_t^{\varphi}\).

As for the numerical realization of the macroscopic equation, we employ a Riemann solver, or more precisely the Harten–Lax–Leer (HLL) Riemann solver [35]. In this method, an approximation for the intercell numerical flux is obtained directly, without the need to solve the local Riemann problems exactly. Therefore this is only an approximate Godunov method. The grid and time step sizes are chosen appropriately to satisfy the CFL condition demanded by the method.

All numerical simulations were implemented in python 2.7.6 with additional scientific computing packages, such as numpy and scipy.

5.1. Comparison with the classical SIR model

Here, we address the question of whether the microscopic model (1) can recover results obtained from the classical SIR model (cf. Section 1), at least in the qualitative sense. Obviously, this would require us to construct an appropriate potential landscape \(\mathcal H\) and interaction force \(\mathcal K\). However, while the microscopic model has multiple functions as parameters, the standard SIR model only has two parameters, namely \(\beta\) and \(\gamma\). For this reason, it is crucial to correctly understand and interpret these parameters accordingly. We further restrict the microscopic model to agents having no mobility \((\sigma=0)\) that are located on an equidistant grid in the domain \(\Omega\). This reduces the model to a deterministic ordinary differential equation, apart from the initial distribution.

As mentioned before, the parameter \(\beta\) in the classical SIR model is known as the transmission rate, which depends on the probability of transmission p and the average number of contacts per agent \(C_0 = C_0(N)\), irregardless of an agent's activity. More specifically, \(\beta = pC_0\). Since the agents are stationary, it is possible to explicitly determine the number of contacts per agent.

Using the indicator function \(\mathbb{1}_{B_R(x)}\), we determine \(C_0\) the number of agents within the vicinity of \(x \in \Omega\) that are maximal radius R away from \(x \in \Omega\) (cf. Fig. 4). On the other hand, if we assume a uniform distribution for the spatial density, i.e., \(\rho_{\text{stat}} \equiv 1/|\Omega|\), then \(C_0\) may be considered as the product of the number density \(N\rho_{\text{stat}}\) and the area of interaction indicated by \(\mathbb{1}_{B_R(x)}\), i.e.,

\[C_0 = C_0(N) = N \int_{\Omega} \mathbb{1}_{B_R(x)}(y) \, d\rho_{\text{stat}}(y) = \frac{N}{|\Omega|} |B_R(x)| = N\pi R^2 / |\Omega|,\]

Fig. 4: Agents on an equidistant grid in \(\Omega\) with possible interaction regions.

which depends explicitly on the number of agents. Considering the equidistant grid in \(\Omega=[0,1]^2\) for any \(N\gg 1\), we rescale the radius as \(R=R_0/\sqrt{N-1}\), where \(R_0>0\) is a fixed constant, in order to keep the number of individual contacts bounded as \(N\to\infty\). More precisely, we have

\[C_0(N) = \frac{N}{N-1} \frac{\pi}{|\Omega|} R_0^2 \longrightarrow c_0 := \frac{\pi}{|\Omega|} R_0^2, \quad \text{as } N \to \infty.\] (13)

There are also other ways of scaling the contact rate (see, for example [21]). However, this consideration is, in fact, the weak coupling scaling mentioned in Remark 1, which allows for the mean-field limit. Indeed, if we set \(\rho^N = \frac{1}{N} \sum_j \delta_{X^j}\) as the empirical measure of locations, then

\[\sum_{j} \mathbb{1}_{B_R(x)}(X^j) = N\rho^N(B_R(x)) \approx \rho^N(B_{R_0}(x)) = \frac{1}{N} \sum_{j} \mathbb{1}_{B_{R_0}(x)}(X^j),\] for a large number of agents \(N\gg 1\). Hence, it makes sense to use the activation function

\[\Phi(x,y) = \mathbb{1}_{B_{R_0}(x)}(y),\] where \(R_0\) is chosen appropriately, depending on \(\beta\).

Fig. 5: Interaction force acting on the class S.

We now work towards identifying the transmission probability p by considering the mean-field equation (2) with the product distribution \(f_t = \rho_{\text{stat}} \otimes g_t\) and \(\psi = \mathbb{1}_{\mathcal{S}}^{\varepsilon}\), \(\chi = c_{\chi} \mathbb{1}_{\mathcal{I}}^{\varepsilon}\), as adopted in (12). The expected effective intensity of interaction between the class of susceptible and infectious agents (cf. Fig. 5) may then be computed as

\[\mathbb{E}\left[\mathcal{F}[f_t](\bar{X}_t,\bar{U}_t)\right] = c_\chi \iint_{S\times S} \mathbb{1}_{B_{R_0}(x)}(y) \mathbb{1}_{\mathcal{S}}^\varepsilon(u) \mathbb{1}_{\mathcal{I}}^\varepsilon(\nu) df_t(y,\nu) df_t(x,u) = c_\chi c_0 S_t I_t,\] where \(S_t = \int \mathbb{1}_{\mathcal{S}}^{\varepsilon} dg_t\) and \(I_t = \int \mathbb{1}_{\mathcal{I}}^{\varepsilon} dg_t\) provides the probabilities in classes \(\mathcal{S}\) and \(\mathcal{I}\), respectively. A direct comparison with the classical SIR model reveals the correspondence \(\beta = c_{\chi} c_0\). However, since \(c_0\) denotes the average number of contacts per agent, we have the relation \(p = c_{\chi}\). Therefore, \(c_{\chi}\) may also be considered as the probability of transmission of a specific disease. Let us summarize the discussion so far. Given \(\beta = pC_0\) from the classical SIR model, we choose \(R_0\) satisfying (13), thereby yielding the interaction term

\[F_N(X^i, U^i, \mathbf{X}, \mathbf{U}) = \frac{p}{N} \sum_{i \neq i} \mathbb{1}_{B_{R_0}(X^i)}(X^j) \mathbb{1}_{\mathcal{S}}^{\varepsilon}(U^i) \mathbb{1}_{\mathcal{I}}^{\varepsilon}(U^j).\]

As for the potential landscape \(\mathcal{H}\), we first note that the classical SIR model does not describe the resistance of an agent towards an infection. Therefore, \(\mathcal{H}'|_{\mathcal{S}} \equiv 0\). On the other hand, if an agent has been infected, i.e., \(U^i \in \mathcal{I} \cup \mathcal{R}\) the interaction term vanishes. Hence, \(\mathcal{H}|_{\mathcal{I} \cup \mathcal{R}}\) should describe the process of recovery. We further assume that \(\mathcal{H}|_{\mathcal{I} \cup \mathcal{R}}\) is linear, with a maximum at \(u = u_*\) and minimum at u = 1 (cf. Fig. 2(b)). Then the evolution of an infected agent is given by

\[\frac{d}{dt}U_t^i = -\mathcal{H}'(U_t^i) = \lambda, \qquad U_0^i = u_*,\] with \(\lambda>0\), where \(-\lambda\) is the slope of \(\mathcal{H}|_{\mathcal{I}\cup\mathcal{R}}\). Solving this equation gives \(U^i_t=u_*+\lambda t\). Recalling the definition of the recovery rate \(\gamma=1/\tau\), where \(\tau>0\) denotes the mean waiting time until an infected individual recovers, we deduce \(U^i_{\tau}=1\). Hence, we obtain the relation

\[\lambda = (1 - u_*)/\tau = \gamma (1 - u_*).\]

Putting all conditions together, we end up with a piecewise linear potential landscape satisfying

\[\mathcal{H}' = \gamma (1 - u_*) \mathbb{1}_{T \sqcup \mathcal{R}}. \tag{14}\]

The piecewise linear potential landscape \(\mathcal{H}_{lin}\) depicted in Fig. 2, for instance, verifies the above requirements and provides a prototype for this comparison throughout this section.

Fig. 6: Compartmental evolution for the microscopic model on an equidistant grid with N=100 agents and M=100 realizations in comparison with classical SIR model. The dashed line represents the classical SIR model.

Fig. 6 shows an acceptable amount of similarity between the two models since the classical SIR model provides solutions that lie within the shaded region of the microscopic model. For the simulation in Fig. 6, we distribute the agents' locations on an equidistant grid, with 90% of the agents having activity uniformly distributed in \(\mathcal{S}\) and 10% of the agents having activity uniformly distributed in \(\mathcal{I}\). Other parameters used are \(C_0=8,\ p=0.5,\ u_*=0,\ \bar{u}=0.3\) and \(\gamma=1.5\).

Remark 5. If we determine \(C_0\) via a standard Gaussian distribution with variance \(\sigma^2 > 0\), instead of using \(\mathbb{1}_{B_R(x)}\), we obtain the formula

\[C_0(N) = N \int_{\Omega} \exp\left(-\frac{|x-y|^2}{2\sigma^2}\right) d\rho_{\rm stat}(y) = 2N\pi\sigma^2/|\Omega|,\] which depends again on the number of agents. To ensure that \(C_0\) remains constant for any \(N \in \mathbb{N}\), we rescale the variance as \(\sigma^2 = \sigma_0^2/N\), which yields \(C_0 = 2\pi\sigma_0^2/|\Omega|\). As opposed to the indicator function, the Gaussian distribution considers every agent as a neighbor. Neighbors that are closer are given more weight than those that are further away. This might be more appropriate whenever considering a domain \(\Omega\) that represents, for example, an enclosed medium sized room.

Fig. 7: Compartmental evolution for the microscopic model on an equidistant grid with N=225 agents, M=100 realizations and a Gaussian-type activation function in comparison with classical SIR model.

Fig.7 depicts the comparison between the microscopic model with a Gaussian-type activation function \(\Phi\), as described above. One observes a slight disparity between the two models, especially in the non-epidemic case, where the number of recovered agents are fewer than the susceptible ones, in contrast to the classical SIR model. Nevertheless, the qualitative behavior of the solutions do coincide to some extend.

5.2. Links between the microscopic and macroscopic models

Another interesting context for the agent-based model is its connection to its macroscopic counterpart. We now investigate this relationship via two ways, namely the direct link, and the mobility link. Unless stated otherwise, we consider a linear potential landscape \(\mathcal H\) satisfying (14), and \(\psi=\mathbb{1}^\varepsilon_{\mathcal S},\ \chi=c_\chi\mathbb{1}^\varepsilon_{\mathcal I}\) for the simulations within this subsection. Furthermore, the agents' locations are initially distributed on an equidistant grid, with 80% of the agents having activity uniformly distributed in \(\mathcal S\) and 20% of the agents having activity uniformly distributed in \(\mathcal I\).

1) Direct link: Looking back at the derivation of the macroscopic model, we first derived the mean-field equation by passing to the limit \(N \to \infty\), and thereafter the spatial activation function \(\Phi\) was removed in the process. Therefore, we will need to look for an appropriate \(\Phi\) for the microscopic model for comparison. In fact, the spatial scaling \(\tilde{x} \sim \varepsilon x\) conducted in Model 3 contracts the bounded domain \(\Omega\) into a spatially concentrated point as \(\varepsilon \to 0\). Hence, the macroscopic model may also be seen as a complete mixture model, where the support of \(\Phi\) is the entire domain, i.e., \(\operatorname{supp}(\Phi) = \Omega\). Consequently, choosing \(\Phi \equiv 1\) results in an agent based model, which is independent of spatial resolution. Since the spatial configuration is obsolete in this case, we may consider any spatial location for the agents.

As seen in Fig. 8, the solution provided by the microscopic model evidently converges to the solution of the macroscopic model as \(N \to \infty\). This verifies on one hand the mean-field limit discussed in Section 3, as well as the choice \(\Phi \equiv 1\). To supplement the validation, we investigate the behavior of the probability distribution corresponding to the microscopic model with N=1600 agents and the macroscopic model, respectively. From this point of view, we recover the complete information concerning the temporal evolution of activity, which provides comprehensive behavior of transitions between the three health classes \(\mathcal{S}\), \(\mathcal{I}\) and \(\mathcal{R}\).

Fig. 8: Compartmental evolution for the microscopic model on an equidistant grid with increasing number of agents in comparison with the evolution generated by the macroscopic model. Other parameters used are C⁰ = 8, c^χ = 0.5, u^∗ = 0, u¯ = 0.3 and γ = 1.5.

Fig. 9: Marginal distribution in the activity space of the microscopic model with N = 1600 agents with parameters as used in Fig. 8.

Fig. 10: Distribution g^t corresponding to the macroscopic model (9) with parameters as used in Fig. 8.

Fig. 11: Convergence of the microscopic model with N = 225 for increasing mobility σ towards the complete mixture model. The dashed line represents the stationary microscopic model, while the dotted line depicts the complete mixture model. Other parameters used are C⁰ = 8, c^χ = 0.5, u^∗ = 0, u¯ = 0.3 and γ = 1.5.

A quick comparison of the figures shown in Fig. 9 and Fig. 10 clearly shows evidence of conformity between the two models, even for the complete distribution sense. Fig. 10 additionally depicts the convergence towards the stationary state g^∞ = δ¹ derived in Section 4.2.

2) Mobility link: In this case, we consider spatial movements of the agents in the form of the standard Wiener process within the bounded domain Ω = [0, 1]² , i.e., we consider the complete microscopic equation (1) with reflective boundary conditions for the Wiener process. The simulations in this part concentrates on the connection between the stationary microscopic model introduced in Section 5.1, which coincides with the classical SIR model, and the complete mixture model described in the previous case. Practically speaking, the presence of mobility 'interpolates' between theses two scenarios, with σ ∈ [0, 1], i.e., the intensity of mobility being the interpolation parameter, as may be seen in Fig. 11.

The results shown in Fig. 11 suggest that increasing mobility increases the mixture of susceptible and infectious agents. This observation is to be expected since the increase in mobility speeds up the rate at which the spatial distribution reaches uniformity, which expresses that the amount of time an agent remains in a particular location is the same amount of time it remains everywhere else. This means that the Wiener process with σ = 1 ensures that every susceptible agent has the same probability to meet an infectious agent as with other susceptible ones, and vice versa. Hence, a microscopic system with a high mobility intensity may just as well be modeled by the macroscopic model which demands less computational effort, since the spatial activation function Φ in the microscopic model has to be recomputed at every time step.

5.3. Microscopic model: A spatially inhomogeneous setting

We finally consider a case the classical SIR model is unable to capture, namely the case where initial distribution of susceptible and infectious agents are no longer uniformly distributed within Ω. To simplify the presentation, we consider an inhomogeneous setting (in initial activity/health status) for the microscopic model with locations distributed equidistantly in Ω.

In Fig. 12, one clearly observes that the distribution of location for infectious agents remains inhomogeneous at time t > 0. Note that the activity of susceptible agents are evenly distributed in activity space. In this particular example, the susceptible agents bordering the upper region of the infectious agents at time t = 0 have an activity close to being immune, i.e., u ≈ −1, and therefore remain susceptible for all times t ≥ 0.

Fig. 13 depicts the compartmental evolution of the microscopic simulation shown in Fig. 12. Plots of this sort cannot be observed when using the classical SIR model, especially in the qualitative behavior of classes I and R between t = 0 and t = 0.5.

5.4. Verification of Theorem 4

Here, we verify the validity of Theorem 4. Therefore, following the assumptions made on the potential landscape and interaction term in Section 4.2, we set

\[\psi = \mathbb{1}_{\mathcal{S}}^{\varepsilon}, \qquad \chi = c_{\chi} \mathbb{1}_{\mathcal{I}}^{\varepsilon}, \qquad \mathcal{H}' = \lambda \mathbb{1}_{\mathcal{S}}^{\varepsilon} - \gamma \mathbb{1}_{\mathcal{I}}^{\varepsilon},\] with S = (−1, 0) and I = (0, 1), for different choices of parameters cχ, λ, γ and initial conditions S0, I0, to obtain the cases R⁰ = ¯ρcχS0/γ − λS0/(γI0) > 1 (epidemic) or R⁰ < 1 (non-epidemic). As in Section 5.2, the agents' locations are initially distributed on an equidistant grid, with 80% of the agents having activity uniformly distributed in S and 20% of the agents having activity uniformly distributed in I. In all cases, we set ρ¯ = 2, to allow for larger values of λ and γ, thereby speeding up the evolution, while keeping c^χ ∈ (0, 1) fixed. We also set λ = γ, which leaves only γ to be varied. In this case, the basic reproduction number simplifies to R⁰ = ¯ρcχS0/γ − S0/I0.

Fig. 14 provides the verification of Theorem 4. In Fig. 14(a), one clearly observes both cases, namely epidemic and non-epidemic, when the basic reproduction number R⁰ is either greater than 1 or less than 1. This figure also affirms the stationary states suggested in Section 4.2. Indeed, for R⁰ > 1, we see that g^t converges towards g^∞ = δ¹ as t → ∞, while for R⁰ < 1, the stationary state is g^∞ = (1 − α)δ⁻¹ + αδ¹ with α ≈ 0.3. Fig. 14(b), on the other hand, provides the evolution of the effective transition Et. As expected the temporal derivative of E^t at t = 0 has the correct sign as indicated by Theorem 4. Essentially, Fig. 14(b) describes the complete behavior of the transition from the class S to class I within an infinitesimal neighborhood of the point u^∗ ∈ J.

Fig. 12: Microscopic model with N = 100 agents on an equidistant grid and a clustered initial distribution of activity. Other parameters used are C⁰ = 8, c^χ = 0.5, u^∗ = 0, u¯ = 0.3 and γ = 1.5.

Fig. 13: Compartmental evolution corresponding to Fig. 12.

6. SUMMARY AND OUTLOOK

In this paper, we successfully developed a model for mathematical epidemiology with spatial resolution, which also allows for a more detailed description of an agent's health status. This clearly paves a way for a more general description of disease transmission, thereby rendering it possible not only to determine the total number of individuals in a certain class of health state, but would also assist in locating the source of a disease. Due to the diversity of the parameters involved in its derivation, the model is able to describe various situations. However, this flexibility becomes also a drawback since the specification of these parameters is not easily accessible and therefore deserves further investigation.

Fig. 14: Evolution of the macroscopic model with \(c_{\chi}=0.3\). The solid lines depict the case \(\Re_0=2\) (epidemic) with \(\gamma=0.08\), while the dashed lines represent the case \(\Re_0=0.8\) (non-epidemic) with \(\gamma=0.1\).

We further provided a recipe for deriving spatial macroscopic models for disease dynamics via passage to the mesoscopic scale. For our specific macroscopic model, we were able to derive a quantity that reflects upon the basic reproduction number \(\mathfrak{R}_0\), which determines the possibility of an outburst. We also showed that stable stationary states exist for the macroscopic equation and that these stationary states are of the form \(g_{\infty} = (1 - \alpha)\delta_{-1} + \alpha\delta_1\), \(\alpha \in [0, 1]\). Unfortunately, the determination of \(\alpha\), depending on the initial distribution \(g_0\) remains an open problem and would therefore require further investigation.

An obvious extension of the current microscopic model would be to incorporate spatial demographic information, as well as spatial interactions among agents. In fact, this was already pointed out in Remark 3. More specifically, one may consider the interacting system

\[dX_t^i = - \left[ \nabla_x V(X_t^i) + \frac{1}{N} \sum\nolimits_{j \neq i} U(X_t^i, X_t^j) \right] dt + \sqrt{2\sigma} \, dW_t^i,\] where \(V: \mathbb{R}^d \to \mathbb{R}\) describes the landscape of an area of a populated region, and \(U: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) is an interaction potential which may be both attracting and repulsive. The convex neighborhood of the local minima of V provide areas that model higher concentration of population, which have a reduced communication between the clusters but still allow for transitions between these clustered populations. This may represent, for example, cities and meeting points. Such potentials may also be used to model special paths on which agents may travel, such as the migration of animals.

Another possible extension is to include different types of agents. Such models may be used to describe vector-based transmitted diseases such as malaria, dengue fever and the Zika virus. In this case, the importance of spatial inhomogeneity is indispensable. The mean-field and macroscopic equations corresponding to such systems are then coupled partial differential equations.

An interesting aspect for modification is the activity set J, which was fixed as J=[-1,1] in this paper. Using a different set J could lead to different dynamics for the activity variable. For instance, we may take \(J=\mathbb{S}^1=\{x\in\mathbb{R}^2\,|\,|x|=1\}\) as the unit circle in \(\mathbb{R}^2\). On this activity set, one can allow for recovered agents to become susceptible again after having been infected, thereby leading to generalization of the well-known classical SIS model. The spatially homogeneous nonlocal macroscopic equation analogous to (9) will then be posed on \(J=\mathbb{S}^1\), or equivalently on [-1,1] with a periodic boundary condition for \(g_t\) on [-1,1], i.e., \(g_t(-1)=g_t(1)\) for all times \(t\geq 0\).

All in all, the basic models introduced in this paper can pave a way to further generalizations that should be numerically investigated and thoroughly analyzed on every scale.

APPENDIX

Without loss of generality, we may suppose \(Z_0^i = \bar{Z}_0^i\), since they are identically distributed. Then, taking the difference of the solutions leads to

\[d(X_t^i-\bar{X}_t^i)=0, \qquad d(U_t^i-\bar{U}_t^i)=-a_N^i\,dt+b_N^i\,dt,\] where the last two terms are given by

\[a_N^i = \mathcal{H}'(U_t^i) - \mathcal{H}'(\bar{U}_t^i),\]
\[b_N^i = \frac{1}{N} \sum_{j \neq i} \mathcal{K}(Z_t^i, Z_t^j) - \int_{\mathcal{S}} \mathcal{K}(\bar{Z}_t^i, z') f_t(dz').\]

The term \(b_N^i\) may be further decomposed to obtain \(b_N^i = c_N^i + d_N^i\), where

\[\begin{split} c_N^i &= \frac{1}{N} \sum\nolimits_{j \neq i} \left[ \mathcal{K}(Z_t^i, Z_t^j) - \mathcal{K}(\bar{Z}_t^i, \bar{Z}_t^j) \right], \\ d_N^i &= \frac{1}{N} \sum\nolimits_{j \neq i} \mathcal{K}(\bar{Z}_t^i, \bar{Z}_t^j) - \int_S \mathcal{K}(\bar{Z}_t^i, z') f(dz'). \end{split}\]

Morever, the sums may be extended to sums over all of j since K is feasible, i.e., K(x, u, y, u) = 0 for any \(x, y \in \Omega\), \(u \in J\). We now investigate the terms separately.

We begin with the term \(a_N^i\) that is easily estimated due to the regularity of \(\mathcal{H}\),

\[-\mathbb{E}[\langle U_t^i - \bar{U}_t^i, a_N^i \rangle] \le c_{\mathcal{H}} \mathbb{E}[|Z_s^i - \bar{Z}_s^i|^2],\] where \(\langle \cdot, \cdot \rangle\) is usual the scalar product on an Euclidean space. As for \(c_N^i\), we simply use the Lipschitz continuity of \(\mathcal{K}\) to obtain

\[\mathbb{E}[\langle U_t^i - \bar{U}_t^i, c_N^i \rangle] \leq c_{\mathcal{K}} \bigg[ \mathbb{E}[|Z_s^i - \bar{Z}_s^i|^2] + \frac{1}{N} \sum\nolimits_{j=1}^N \mathbb{E}[|Z_s^j - \bar{Z}_s^j|^2] \bigg].\]

To estimate \(d_N^i\), we first define

\[\kappa_t^i(\bar{Z}_t^j) := \mathcal{K}(\bar{Z}_t^i, \bar{Z}_t^j) - \mathbb{E}[\mathcal{K}(\bar{Z}_t^i, \bar{Z}_t^j)] = \mathcal{K}(\bar{Z}_t^i, \bar{Z}_t^j) - \int_S \mathcal{K}(\bar{Z}_t^i, z') f_t(dz'),\]

Clearly \(\mathbb{E}[\kappa_t^i(\bar{Z}_t^j)|\bar{Z}_t^i]=0\) for any \(j\neq i\). Furthermore, we have

\[\mathbb{E}[\kappa_t^i(\bar{Z}_t^j)\kappa_t^i(\bar{Z}_t^k)] = \mathbb{E}[\mathbb{E}[\kappa_t^i(\bar{Z}_t^j)\kappa_t^i(\bar{Z}_t^k)|\bar{Z}_t^i]] = \mathbb{E}[\mathbb{E}[\kappa_t^i(\bar{Z}_t^j)|\bar{Z}_t^i]\mathbb{E}[\kappa_t^i(\bar{Z}_t^k)|\bar{Z}_t^i]] = 0,\]

since the processes \(\bar{Z}_t^j\) and \(\bar{Z}_t^k\) are independent for \(j \neq k \neq i\). Consequently

\[\mathbb{E}\left[\left|\frac{1}{N}\sum_{j\neq i}\kappa_t^i(\bar{Z}_t^j)\right|^2\right] = \frac{1}{N^2}\sum_{j,k\neq i}\mathbb{E}\left[\kappa_t^i(\bar{Z}_t^j)\kappa_t^i(\bar{Z}_t^k)\right] = \frac{N-1}{N^2}\mathbb{E}\left[\kappa_t^i(\bar{Z}_t^j)^2\right] \\ \leq \frac{N-1}{N^2}\iint_{S\times S}|\mathcal{K}(z,z')|^2f_t(dz')f_t(dz).\]

Therefore, we obtain, by Young's inequality, the estimate

\[\begin{split} \mathbb{E}[\langle Z_t^i - \bar{Z}_t^i, d_N^i \rangle] &\leq \frac{1}{2} \mathbb{E}[|Z_t^i - \bar{Z}_t^i|^2] + \frac{1}{2} \mathbb{E}[|d_N^i|^2] \\ &\leq \frac{1}{2} \mathbb{E}[|Z_t^i - \bar{Z}_t^i|^2] + \frac{N-1}{2N^2} \iint_{S \times S} |\mathcal{K}(z, z')|^2 f_t(dz') f_t(dz) \\ &\leq \frac{1}{2} \mathbb{E}[|Z_t^i - \bar{Z}_t^i|^2] + \frac{c_0}{N}. \end{split}\]

Now, set \(Y_t^i = \mathbb{E}[|Z_t^i - \bar{Z}_t^i|^2]\). Then, by Itô's calculus and the estimates above, we obtain

\[\frac{d}{dt}Y_t^i = 2 \operatorname{\mathbb{E}}\left[ \langle Z_t^i - \bar{Z}_t^i, a_N^i + c_N^i + d_N^i \rangle \right] \leq c \left[ Y_t^i + \frac{1}{N} \sum\nolimits_{j=1}^N Y_t^j + \frac{1}{N} \right].\]

Averaging over 1 ≤ i ≤ N gives

\[\frac{d}{dt}Y_t^{(N)} := \frac{d}{dt}\frac{1}{N}\sum\nolimits_{i=1}^{N}Y_t^i \le \tilde{c}\left[Y_t^{(N)} + \frac{1}{N}\right].\]

An application of the Gronwall inequality yields

\[Y_t^N \le \frac{\tilde{c}}{N} t e^{\tilde{c}t}.\]

Substituting this into the inequality for Y i t and using Gronwall's inequality again yields

\[\sup_{t \in [0,T]} Y_t^i \le \frac{\tilde{c}}{N} T e^{2\tilde{c}T},\] which is precisely the required estimate for any 1 ≤ i ≤ N ∈ N.