From interacting agents to density-based modeling with stochastic PDEs

TL;DR

This paper develops a stochastic PDE-based density model from agent-based models to efficiently simulate large populations, with applications demonstrated through innovation spreading examples.

Contribution

It introduces a novel reduction from agent-based models to stochastic PDEs, enabling scalable simulation of large systems with regularization via finite element discretization.

Findings

SPDE model effectively approximates agent-based models for large populations

Finite element discretization improves simulation efficiency and regularizes the SPDE

Illustrative examples demonstrate the model's applicability to innovation spreading

Abstract

Many real-world processes can naturally be modeled as systems of interacting agents. However, the long-term simulation of such agent-based models is often intractable when the system becomes too large. In this paper, starting from a stochastic spatio-temporal agent-based model (ABM), we present a reduced model in terms of stochastic PDEs that describes the evolution of agent number densities for large populations. We discuss the algorithmic details of both approaches; regarding the SPDE model, we apply Finite Element discretization in space which not only ensures efficient simulation but also serves as a regularization of the SPDE. Illustrative examples for the spreading of an innovation among agents are given and used for comparing ABM and SPDE models.

Figures

Equations

$$d{X}_i(t)=-\nabla V(X_i(t))dt+\sigma d{B}_i(t),$$

$$d{X}_i(t)=-\left(\nabla V(X_i(t))+\nabla U_i(X(t))\right)dt+\sigma d{B}_i(t).$$

$$U_i:D^N\subseteq {\mathbb{R}}^{d\times N}\to \mathbb{R},X(t)\mapsto \sum _{j\ne i}u\left(\Vert X_i(t)-X_j(t)\Vert \right),$$

$$R_r:T_s+T_{s^{\prime \prime }}\to T_{s^{\prime }}+T_{s^{\prime \prime }}\text{with}s,s^{\prime },s^{\prime \prime }\in \{1,\dots ,N_T\}$$

$${\lambda }_i^r(t)={\lambda }_i^r(\mathcal{A}(t),Y(t))={\gamma }_r\sum _{j=1}^N{\mathcal{A}}_{ij}(t)\text{mathbbm}{1}_{\{s^{\prime \prime }\}}(Y_j(t))\text{mathbbm}{1}_{\{s\}}(Y_i(t)),$$

$$Y_i(t)=Y_i(0)+\sum _{r=1}^{N_R}{\mathcal{P}}_i^r\left({\int }_0^t{\lambda }_i^r(t^{\prime })d{t}^{\prime }\right)v_r,$$

$$X_i(t)$$

$$Y_i(t)$$

$${\mathcal{R}}_s:=\{r\in \{1,...,N_R\}:R_r:T_s+T_{s^{\prime \prime }}\to T_{s^{\prime }}+T_{s^{\prime \prime }}\}$$

$$X_i(t_{k+1})=X_i(t_k)-\nabla V(X_i(t_k))\Delta t+\sigma \sqrt{\Delta t}{\zeta }_{i,k}$$

$${\Lambda }_i(t_k):=\sum _{r\in {\mathcal{R}}_s}{\lambda }_i^r(t_k)$$

$$1-\exp (-{\Lambda }_i(t_k)\Delta t).$$

$${\lambda }_i^r(t_k)/{\Lambda }_i(t_k)$$

$${\rho }_s:D\times [0,T]\mapsto {\mathbb{R}}_{\ge 0}$$

$$\frac{\partial {\rho }_s(x,t)}{\partial t}=\mathcal{D}{\rho }_s(x,t)+\mathcal{I}{\rho }_s(x,t)$$

$$\mathcal{D}{\rho }_s(x,t):=\frac{{\sigma }^2}{2}\Delta {\rho }_s(x,t)+\nabla \cdot (\nabla V(x){\rho }_s(x,t))+\sigma \nabla \cdot \left(\sqrt{{\rho }_s(x,t)}{Z}_s^{\mathcal{D}}(x,t)\right)$$

$$\mathbb{E}\left(Z_{s,j}^{\mathcal{D}}(x,t)Z_{s^{\prime },j^{\prime }}^{\mathcal{D}}(x^{\prime },t^{\prime })\right)={\delta }_{j{j}^{\prime }}{\delta }_{s{s}^{\prime }}\delta (x-x^{\prime })\delta (t-t^{\prime }),$$

$${\rho }^i(x,t)=\delta (x-X_i(t)),$$

Full text

Keywords Agent-based modeling, Model reduction, Dean-Kawasaki model, SPDEs, Finite Element method

1 Introduction

Modeling real-world dynamics is of great importance for the understanding, prediction and manipulation of the dynamical process of interest. Mathematically, one is therefore interested in the model analysis, inference of its parameters, as well as the accurate simulation and control of the dynamics. The modeling of realistic processes is nowadays often based on formulations in terms of discrete and autonomous entities (e.g., humans, institutions, companies), so-called agents, that move in a given environment and interact with each other according to a set of rules. Agent-based models (ABMs) are very flexible in their description, often they just consist of a set of behavioural rules for agents encoded in computer programs [bonabeau2002agent, macy2002factors, galan2009errors]. Due to their versatility and the ease of incorporating data and assumptions into the model, ABMs are very attractive for modelers of various complex applications and disciplines (e.g., social sciences [bankes2002agent, macy2002factors], ecology [grimm2005pattern, grimm2013individual], archaeology [griffith2010hominids, graham2006networks]). However, the different existing formulations of ABMs appear to be quite inconsistent and lack rigorous theoretical foundation, which motivates to develop standardized formats for these models [grimm2010odd, bankes2002agent, gilbert2002platforms].

In [ourepjpaper, conrad2018mathematical], a mathematical formulation of an ABM is presented for systems of spatially distributed agents that move randomly in space and interact whenever they are close-by. The local interactions among agents trigger them to change their type, e.g., their opinion, or their infection status. The model is formulated as a system of equations coupling the Markovian diffusion dynamics for the spatial agent movement with Markov jump processes for the type changes of each agent. Agent-based models of this class can for instance be found as models for infection spreading [brockmann2013hidden, pastor2015epidemic], innovation spreading [ourepjpaper, conrad2018mathematical] and chemical reactions [doi1976stochastic], but also other classes of ABMs are commonly formulated as Markov processes [izquierdo2009techniques].

In general, these ABMs cannot be solved analytically, but require numerical methods for simulating realizations of the modeled dynamics. In [ourepjpaper], a joint algorithm for the simulation of the spatial diffusion and the event-based interaction process is presented, building on the Euler-Maruyama scheme [higham2001algorithmic] and the Temporal Gillespie algorithm [vestergaard2015temporal]. As the model formulation is stochastic, many Monte-Carlo (MC) simulations are required in order to make adequate predictions for observables of the system. However, for most real-world dynamics, such MC simulations are intractable due to an explosion of the costs for increasing agent numbers.

In Figure 1, we plot one realization of a MC simulation for innovation spreading in ancient times [ourepjpaper, conrad2018mathematical]. The underlying ABM model was built to study the spreading of the wool-bearing sheep in the Near East and Southeast Europe between $6000$ and $4000$ BC. This stochastic model considers $N=4000$ agents that diffuse in a complex, data-driven suitability landscape with non-trivial boundary conditions. Due to the high computational costs of simulating this model and the need for running many such MC simulations, developing new approaches for model reductions with a small approximation error is essential.

For spatial systems of many interacting agents, one can reduce the model complexity by replacing the micro-scale model of individual agents by a meso-scale model of stochastic agent densities, leading to a system of coupled stochastic partial differential equations (SPDEs) [dean1996langevin, kawasaki1973new, bhattacharjee2015fluctuating, kim2017stochastic]. Classically, in the limit of infinitely many agents, stochastic effects become negligible and the dynamics can be described by a macro-scale model of deterministic reaction-diffusion PDEs [pearson1993complex, fisher1937wave]. In a real-world system, however, the number of agents is typically finite, such that the dynamics are still intrinsically random and stochastic modeling approaches are more appropriate in order to reflect the uncertainty and fully cover the emerging phenomena of the underlying ABM.

In this paper, after reviewing and slightly generalizing the stochastic ABM from [ourepjpaper, conrad2018mathematical], we will formulate a system of SPDEs as a meso-scale approximation to the ABM, see [helfmann2018stochastic]. The model is an extension of the Dean-Kawasaki model [dean1996langevin, kawasaki1973new] that describes the transport in space of the stochastic agent densities and includes interactions between the agent densities for the different types [kim2017stochastic, bhattacharjee2015fluctuating]. The derived distribution-valued equations belong to the class of Dean-Kawasaki type problems, for which well-posedness is still unresolved due to several mathematical difficulties [cornalba2018regularised, lehmann2018dean, fehrmanwell].

Since in the end we are interested in the discretization and simulation of the model, we will, in order to deal with these issues, interpret the system of SPDEs in its weak form and consider “a regularization by discretization”. Henceforth, we propose an algorithm for the efficient sampling of trajectories of the coupled SPDE system by first discretizing in space using the Finite Element (FE) method, thereby approximating the unknown distributions by piece-wise polynomial functions. The resulting system of SDEs is well-posed and can be further discretized in time using the Euler-Maruyama method.

As a meso-scale approach, simulation of the system of SPDEs is much faster than of the corresponding ABM, since the number of coupled equations is drastically reduced from scaling with the number of agents to scaling with the number of different agent types. Further, the Finite Element method can easily deal with complicated domain boundaries [de2015finite], which is crucial for many agent-based models and real-world situations (see Figure 1). Other approaches, e.g., using Finite Volume schemes have been intensively studied [bhattacharjee2015fluctuating, kim2017stochastic] and benefit from the property of agent number conservation.

Last, we will study the ABM and the reduced SPDE model numerically on a toy example by comparing the computational effort and investigating the approximation quality of the reduced SPDE model with respect to the ABM. The focus of the present paper is on the mapping from agent-based models to SPDE models and their approximation quality and not on a numerical analysis of the approximation, still we will also shortly reflect on the consistency of numerical solutions to the SPDE.

2 Agent-based Model

Agent-based models are micro-scale computational models describing the system of interest in terms of discrete entities, called agents. An agent represents an autonomous entity such as a chemical particle, a person, a group of people or an organization. The individual behaviour of each agent is traced in space and time, including its interactions with other agents and the environment. The hope is that from the interplay of the local behaviour of agents on the microscopic scale, global patterns emerge on a larger system scale  [helbing2012agent, macy2002factors, grimm2005pattern]. Emergence is likewise often stated as “the whole is more than the sum of its parts” due to the interactions among the parts (i.e., the agents) resulting in global structures [bonabeau2002agent].

In this section, we will review and slightly generalize the agent-based model from  [ourepjpaper, conrad2018mathematical], which is formulated as a system of stochastic processes describing the position dynamics and the interaction rules for agents. The position dynamics of agents are described by Brownian motion with a drift term pushing them towards suitable parts of their local surroundings. This is a useful description in the case that agents move randomly and slightly jittery in space whilst following some (underlying) energy landscape, e.g., when describing the mobility of humans in the ancient world, the movement of animals, or the diffusion of chemical particles [doi1976stochastic].

Further, agents can interact according to a set of predefined rules whenever they are close in space. These interactions cause agents to change their type influenced by the state of other near-by agents (e.g., their infection status [kermack1927contribution, pastor2015epidemic], their current opinion [holme2006nonequilibrium], whether they use an innovation [conrad2018mathematical, ourepjpaper], whether they have a certain information). The type changes are commonly written as $T_{1}+T_{3}\rightarrow T_{2}+T_{3}$ (inspired by the notation for chemical reactions), where $T_{1},T_{2},T_{3}$ denote different agent types.

Even though the model formulation is very simple, it generates a vast range of different and complex outcomes on a larger system scale. The presented ABM is based on the Doi model  [doi1976stochastic] for reacting and locally diffusing chemical particles. It has similarities to the SIR model for infection spreading [kermack1927contribution, pastor2015epidemic] and to Brownian agents [schweitzer2002modelling] and builds on the ABM for innovation spreading in the ancient world as introduced in [conrad2018mathematical, ourepjpaper] and Figure 1.

We will give a detailed description of the agent-based model in Section 2.1, and explain a method for its simulation in Section 2.2.

2.1 Model Formulation

We consider a system of $N$ agents, where each agent $i=1,\dots,N$ is characterized by its type $Y_{i}(t)$ and its position $X_{i}(t)$ at time $t\in[0,T]$. The agents’ positions can be continuous values in a given domain $D\subseteq\mathbb{R}^{d}$, whereas the type of an agent is a discrete feature denoted by values in $\{1,\dots,N_{T}\}$. Thus the state of the $i^{\text{th}}$ agent at time $t$ is given by $(X_{i}(t),\ Y_{i}(t))\in D\times\{1,\dots,N_{T}\}.$ We write the state of the whole system of $N$ agents as $(X(t),Y(t))=(X_{i}(t),\ Y_{i}(t))_{i=1}^{N}.$

In the model, we are following every agent $i=1,\dots,N$ individually and track the evolution in time of its position state and type. The position dynamics and the interactions between agents leading to type changes are described in the following sections.

2.1.1 Position Dynamics of Agents

We assume that agents are able to move and change their position in the domain $D\subseteq\mathbb{R}^{d}$. Thereby agents are taking their local surroundings into account, in such a way, that they are attracted to near-by regions that are suitable for them and refrain from unsuitable parts of the domain.

We can straightforwardly model this by letting the agents follow the gradient of a potential landscape, the so-called suitability landscape. The suitability landscape $V$ indicates the attractivity of the environment and gives an incentive for agents to prefer or avoid certain near-by parts of the domain. Valleys of the suitability landscape correspond to attractive regions and peaks and divides correspond to unsuitable areas that are moreover difficult to overcome. The suitability landscape can be constructed on the basis of data and expert knowledge, see e.g., [ourepjpaper, conrad2018mathematical]. Additionally, we include randomness in the agents’ motion to account for other unknown incentives for positional changes and to allow agents to be exploratory or make mistakes in their evaluation of the environment.

Thus, the change of the position $X_{i}(t)\in D\subseteq\mathbb{R}^{d}$ of every agent $i=1,\dots,N$ is governed by the Itô diffusion process

[TABLE]

with gradient operator $\nabla=\left(\frac{\partial}{\partial x_{1}},\dots,\frac{\partial}{\partial x_{d}}\right)^{T},$ suitability landscape $V:D\subseteq\mathbb{R}^{d}\rightarrow\mathbb{R}$, diffusion constant111Note that within physical sciences the term diffusion constant or diffusion coefficient typically refers to $\frac{\sigma^{2}}{2}$, however, in the theoretical context of stochastic differential equations these terms are also used for $\sigma$ itself [oksendal2013stochastic, evans2012introduction]. $\sigma\in\mathbb{R}$ and $B_{i}(t)$ denoting independent standard Brownian motions in $\mathbb{R}^{d}$. We impose reflecting boundary conditions in the case where $D$ is bounded and thereby ensure that the agents’ positions are within $D$ for all times.

Remark*.*

The position dynamics of agents in many modeling scenarios are interdependent such that agents tend to group together in space and form clusters, while also keeping some distance from each other in order to avoid spatial overlap and crowding [conrad2018mathematical, ourepjpaper]. To account for this, another drift term has to be added to the position SDE (1), s.t. the dynamics are governed by

[TABLE]

Thereby each agent $i$ experiences an additional force derived from the attraction-repulsion potential $U_{i}$

[TABLE]

where we sum over all pair-wise attraction-repulsion potentials $u:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}$ between agent $i$ and $j\neq i$. The pair-wise potentials $u$ are inspired by interatomic potentials from Physics (e.g., Lennard-Jones potential, Buckingham potential). Attraction between pairs of agents occurs whenever agents at long distances are driven towards another, and repulsion appears when agents are forced apart at short distances. Agents are thus searching for an optimal balance between forming clusters of agents on the one hand and distributing in space on the other hand.

2.1.2 Interaction Rules

Agents can change their type according to a set of $N_{R}$ interaction rules $\{R_{r}\}$, $r=1,\dots,N_{R}$. The type changes happen at a certain rate and whenever an agent is in proximity of specific other agents that can influence the agents’ type. We consider rules $\{R_{r}\}$ that can be written as the type change222In principle, the interaction rules could be of a more complicated form, e.g., by including the death and birth of agents such as $T_{s}\rightarrow\emptyset$ and $\emptyset\rightarrow T_{s}$, by allowing both agents to change their type $T_{s}+T_{s^{\prime\prime}}\rightarrow T_{s^{\prime}}+T_{s^{\prime\prime\prime}}$ or by letting a larger number of agent types take part in an interaction. A straightforward mapping of these more complicated interaction rules is also possible for the SPDE formulation that will be introduced in Section 3.1, we keep these simple interaction rules only for notational simplicity.

[TABLE]

happening at the fixed influence rate $\gamma_{r}$ and triggered by a close-by agent of type $T_{s^{\prime\prime}}$.

For example when modeling the spreading of an innovation among humans  [ourepjpaper, conrad2018mathematical], we can assume that every agent takes one of the two discrete innovation states: $T_{1}$ for a non-adopter or $T_{2}$ for an adopter of the innovation. By defining one simple interaction rule $R_{1}$: $T_{1}+T_{2}\rightarrow 2\ T_{2}$, we can model that adopters pass on the innovation to non-adopters at a fixed rate whenever they are in contact.

Pairs of agents are in contact if they are within a distance $d_{\text{int}}$ of each other. Given the changing positions of agents (1), we construct a time-evolving network between agents (represented by the set of nodes) that are in contact (given by the edges) at time $t$. The network is fully determined by a time-dependent adjacency matrix $\mathcal{A}(t)=(\mathcal{A}_{ij}(t))_{i,j=1}^{N}$ with entries $\mathcal{A}_{ij}(t)=1$ if $i\neq j,\ \|X_{i}(t)-X_{j}(t)\|\leq d_{\text{int}}$, and 0 else.

As a next step, we are interested in describing the type change process $\{Y_{i}(t)\}_{t\in[0,T]}$ of each agent $i$. If agent $i$ at time $t$ is of type $T_{s}$, $s\in\{1,\dots,N_{T}\}$, we denote this by $Y_{i}(t)=s$. Type changes for agent $i$ are modeled as Markov jump processes on $\{1,\dots,N_{T}\}$ with time-dependent transition rates. The transition rates are changing in time since they depend both on the proximity of other agents and their types. The transition rate function $\lambda_{i}^{r}(t)$ gives the accumulated rate for agent $i$ of type $T_{s}$ to change its type according to interaction rule $R_{r}$ and is proportional to the constant influence rate $\gamma_{r}$ and to the number of neighbors of agent $i$ that trigger interaction $R_{r}$, i.e., the number of agents of type $T_{s^{\prime\prime}}$ in rule (4).

Consequently, we write the transition rate function for agent $i$ at time $t$ following interaction rule $R_{r}:\ T_{s}+T_{s^{\prime\prime}}\rightarrow T_{s^{\prime}}+T_{s^{\prime\prime}}$ as

[TABLE]

where $\mathbbm{1}_{B}$ is the indicator function defined as $\mathbbm{1}_{B}(x)=1$ if $x\in B$ and 0 else. Finally, the type change process $\{Y_{i}(t)\}_{t\in[0,T]}$ of agent $i$ can be expressed as

[TABLE]

where $Y_{i}(0)$ denotes the initial type of agent $i$, $\mathcal{P}_{i}^{r}$ denote i.i.d. unit-rate Poisson processes and the type change vector is denoted by $v=(v_{r})_{r=1}^{N_{R}}$ ($v_{r}=s^{\prime}-s$ for $R_{r}$ as given in (4) above).

2.1.3 Formulation of the Agent System Dynamics

Putting together Equations (1) and (6), we can describe the evolution in time of the agent states by the following coupled equations

[TABLE]

for agents $i=1,\dots,N$.

Since we cannot solve the coupled equations (2.1.3) analytically, we will in the following explain how to accurately discretize and efficiently simulate trajectories of the dynamics.

2.2 Simulation

For each agent $i$, the type change process depends on the positions and on the types of all other agents via the time-dependent transition rate functions $\lambda_{i}^{r}(t)$. If we additionally include attraction-repulsion forces between agents, then further the motion of all agents is intrinsically intertwined. For the design of the simulation algorithm we have to take this into account and therefore propose to simulate both the position and type changes (2.1.3) in parallel.

The most straightforward approach is to discretize time $[0,T]$ with sufficiently small $\Delta t$ time steps. The Euler-Maruyama method can be employed to accurately and efficiently discretize the position SDE [kloeden, ourepjpaper], whereas when discretizing the Markov jump processes we replace rates by probabilities in each time step [conrad2018mathematical, ourepjpaper]. In contrast to time-continuous simulation schemes for the interactions, like the standard Gillespie algorithm, such a time discretization can evoke interaction conflicts in each integration step, namely when two type change events are chosen to occur for the same agent. In order to circumvent such conflicts, we define for each type $T_{s}$ the index set

[TABLE]

of rules that change this type, and stipulate that each agent of type $T_{s}$ can undergo maximally one event $R_{r}$, $r\in\mathcal{R}_{s}$, per time step. The resulting simulation approach can be summarized as follows.

For each time step $t_{k}=k\Delta t$, $k=0,\dots,K-1$:

For each agent $i=1,\dots,N$: the positions are advanced according to

[TABLE]

with i.i.d. $\zeta_{i,k}\sim\mathcal{N}(0,1)$.

• 2.

For each $i=1,\dots,N$, set $s:=Y_{i}(t_{k})$ (the type of agent $i$ at time $t_{k}$) and calculate

[TABLE]

with $\mathcal{R}_{s}$ defined in (8). Then, a type change of agent $i$ occurs with probability

[TABLE]

In case the agent changes its type, decide subsequently by which of the rules, choosing $R_{r}$ with probability

[TABLE]

for each $r\in\mathcal{R}_{s}$.

A few remarks are in order: Note that this procedure does not restrict the total number of interaction events per agent and time step to one. An agent may still undergo several interaction events in parallel as long as maximally one of these events changes its type (i.e., for an agent of type $T_{s^{\prime\prime}}$ several events of the form (4) are possible). Further, since the underlying interaction dynamics are time-continuous, a more accurate simulation [fennell2016limitations] of the type changes would simulate them continuously in time using a Gillespie algorithm [vestergaard2015temporal], for the coupled simulation approach in this case we refer the reader to [ourepjpaper].

Regarding the computational cost, we note that the algorithm scales more than linearly with the number of agents $N$, since in each time step, we need to compute the contact network $\mathcal{A}(t_{k})$, i.e., the pair-wise distances between agents. The brute-force approach of computing all pair-wise distances scales like $\mathcal{O}(N^{2})$ for $N$ agents, whereas clever algorithms can lower the cost to $\mathcal{O}(N^{5/3})$  [verlet1967computer]. These distances are also required when including attraction-repulsion forces, since the pair-wise attraction-repulsion potential $u\left(\lVert X_{i}(t)-X_{j}(t)\rVert\right)$ depends on the pair-wise distances between agents.

The proposed simulation algorithm produces a single trajectory of the joint stochastic process $\{X(t),Y(t)\}_{t\in[0,T]}$. In order to obtain reliable statistics of the dynamics and to estimate the process’ distribution and its moments, several thousand repeated simulations are required. In case of a system with large numbers of agents this causes problems since the simulation cost scales badly with increasing agent numbers $N$. Parameter studies are numerically not tractable in this case. Many real-world systems actually contain very large numbers of agents, which motivates to consider a less complex modeling approach in terms of stochastic PDEs, and thereby offers a cheap method for simulating trajectories.

3 Model Reduction to a System of Stochastic PDEs

Instead of studying the evolution of every individual agent, in the following we study the transport in space and the interactions between agents in terms of number densities for each type, such that the model complexity and simulation costs are drastically reduced. This reduced model, given in the form of a system of stochastic PDEs, is a combination of the model for position changes proposed by Dean and Kawasaki [dean1996langevin, kawasaki1973new] (also termed Dean-Kawasaki model), and added terms for the interactions between agent densities [bhattacharjee2015fluctuating, kim2017stochastic]. It approximates the ABM for systems of many, but finitely many agents. Stochasticity still emerges from the systems’ inherent randomness due to the finite number of agents. But agents become indistinguishable among their type and we lose the individual agent labels, see Figure 2.

Applications of similar SPDE models include chemical pattern formation [bhattacharjee2015fluctuating, kim2017stochastic] and models for bacterial populations [grafke2017spatiotemporal], but to our knowledge they have not yet been applied for modeling systems of humans.

We will start by proposing and explaining the features of the reduced SPDE model in Section 3.1. In Section LABEL:sec:discr, we will tackle its efficient numerical discretization by means of the Finite Element method.

3.1 Model Formulation

We are considering a model describing the stochastic evolution of agent number densities (or number concentrations) for each agent type $s=1,\dots,N_{T}$. The (stochastic) agent number density

[TABLE]

is defined on the domain $D\subseteq\mathbb{R}^{d}$, the time interval $[0,T]$ and with probability space $(\Omega,\mathcal{F},\mathbb{P})$333Actually the number density should be $\rho_{s}(\omega,x,t)$ but we do not write the explicit dependence on $\omega\in\Omega$. The agent number density is random since it solves a stochastic PDE.. Integrating the number density $\rho_{s}$ over the domain yields the number of agents of type $T_{s}$, which we denote by $N_{s}$.

The densities $(\rho_{1}(x,t),\dots,\rho_{N_{T}}(x,t))$ evolve due to diffusion and drift in the suitability landscape and because of the set of interaction rules. The temporal changes of $\rho_{s}(x,t)$ for $s=1,\dots,N_{T}$ are given by the following stochastic partial differential equation (SPDE) interpreted in the Itô sense

[TABLE]

with stochastic diffusion operator $\mathcal{D}$ (see Section 3.1.1) and stochastic interaction operator $\mathcal{I}$ (Section LABEL:sec:int_op). For a fixed sample $\omega\in\Omega$, the agent density is a realization of a stochastic process solving the SPDE.

3.1.1 Diffusion of the Agent Densities

The diffusion operator in Equation (9) is given by [dean1996langevin, kawasaki1973new, kim2017stochastic, bhattacharjee2015fluctuating]

[TABLE]

with diffusion constant $\sigma\in\mathbb{R}$, Laplace operator ${\Delta=\sum_{l=1}^{d}\frac{\partial^{2}}{\partial x_{l}^{2}}}$ and suitability landscape ${V:D\subseteq\mathbb{R}^{d}\mapsto\mathbb{R}}$. $Z_{s}^{\mathcal{D}}(x,t)=\left(Z_{s,1}^{\mathcal{D}}(x,t),\dots,Z_{s,d}^{\mathcal{D}}(x,t)\right)$ denotes a $d$-dimensional vector of space-time white noise (STWN) for the diffusion, i.e., a mean zero process that is uncorrelated in space and time

[TABLE]

where $\delta_{ij}$ denotes the Kronecker Delta and $\delta(x-y)$ denotes the Dirac Delta distribution.

The diffusive part of the SPDE evolves a number density of many agents and is responsible for the diffusive transport in space with drift in the suitability landscape $V(x)$. Its deterministic part $\partial_{t}\rho_{s}=\frac{\sigma^{2}}{2}\Delta\rho_{s}+\nabla\cdot(\nabla V\rho_{s})$ reminds of the Fokker-Planck equation that describes the evolution of a density of infinitely many non-interacting particles that are diffusing with drift in a potential. Here though, we consider a large but finite number of diffusing particles, as such the system is still intrinsically random and the derivation leads to an additional noise term fluctuating around zero. The derivation [dean1996langevin, kawasaki1973new] of $\frac{\partial\rho_{s}}{\partial t}=\mathcal{D}\rho_{s}$ from the position SDE (1) is based on the following idea: We define the empirical agent density for a single agent $i$ as

[TABLE]

where $\delta(x-X_{i}(t))$ denotes a Dirac Delta distribution placed at position $X_{i}(t)$. Then, using Itô’s Formula, the SDE (1) for agent $i$ is transformed into an SPDE describing the temporal evolution of $\rho^{i}(x,t)$. We sum the SPDEs over $i=1,\dots,N_{s}$ to get an equation for the density of all agents of type $T_{s}$, i.e., for $\rho_{s}(x,t)=\sum_{i=1}^{N_{s}}\rho^{i}(x,t).$ Aiming at a closed-form equation for $\rho_{s}(x,t)$ that does not depend on $\rho^{i}(x,t)$ explicitly, Dean [dean1996langevin] in his formal computations approximates the obtained conditional Gaussian noise term by a Gaussian noise of the same mean and covariance that just depends on $\rho_{s}(x,t)$.

This multiplicative noise term $\sigma\nabla\cdot(\sqrt{\rho_{s}}Z_{s}^{\mathcal{D}})$ is non-linear and comes with several mathematical problems, making the question of existence and uniqueness of the solution still an open problem [cornalba2018regularised, lehmann2018dean]. First, the noise term is given in a divergence form. In order to interpret this term, we will consider the weak formulation and use partial integration to transfer the divergence operator to the test functions. Second, the meaning of the term $\sqrt{\rho_{s}}$ remains undefined. Namely, the unknown $\rho_{s}$ is defined as a sum of Dirac Delta distributions, hence the meaning of its square root is unclear. We will deal with this problem by approximating the distribution $\rho_{s}$ by a function from a Finite Element space, for which the square root is well defined.

Different notions of solutions to (10) that would be natural to consider exist, such as a martingale solution [lehmann2018dean] or a path-wise kinetic solution [fehrmanwell]. In order to apply the Finite Element method and to reformulate the divergence of the noise term, in our setting it is more appropriate to consider the weak formulation (in the PDE sense) that will be derived in LABEL:sec:weak.