Computational framework to capture the spatiotemporal density of cells with a cumulative environmental coupling

TL;DR

This paper introduces a new computational framework that models cell density considering environmental factors, providing a more efficient alternative to stochastic agent-based simulations for studying cellular processes.

Contribution

The work develops a continuum PDE model coupled with environmental variables, including stability analysis and a numerical method using fundamental solutions.

Findings

Validated the PDE model against agent-based simulations.

Demonstrated the model's ability to capture cell death due to toxic exposure.

Provided a stable and well-posed mathematical formulation.

Abstract

Stochastic agent-based models can account for millions of cells with spatiotemporal movement that can be a function of different factors. However, these simulations can be computationally expensive. In this work, we develop a novel computational framework to describe and simulate stochastic cellular processes that are coupled to the environment. Specifically, through upscaling, we derive a continuum governing equation that considers the cell density as a function of time, space, and a cumulative variable that is coupled to the environmental conditions. For this new governing equation, we consider the stability through an energy analysis, as well as proving uniqueness and well-posedness. To solve the governing equations in free-space, we propose a numerical method using fundamental solutions. As an application, we study a cell moving in an infinite domain that contains a toxic chemical, where a cumulative exposure above a critical value results in cell death. We illustrate the validity of this new modeling framework and associated numerical methods by comparing the density of live cells to results from the corresponding agent-based model.

Figures

Equations

$$f{\ast }^{m}g\equiv \underset{m\text{ times}}{\underbrace{f\ast f\ast \cdots \ast f}}\ast g.$$

$$\begin{array}{rl}U(x,t+& △t,\xi +\hat{\beta }(x))=\ell (x+△x)U(x+△x,t,\xi )\\ & +r(x-△x)U(x-△x,t,\xi )+[1-r(x)-\ell (x)]U(x,t,\xi ),\end{array}$$

$$U(x,t+△t,\xi +\hat{\beta }(x))$$

$$+\mathcal{O}(△t^2,\hat{\beta }(x{)}^2),$$

$$U(x\pm △x,t,\xi )$$

$$+\mathcal{O}(△x^3),$$

$$q(x+△x)$$

$$U+△t\frac{\partial U}{\partial t}+\hat{\beta }(x)$$

$$+[rU]-△x\frac{\partial (rU)}{\partial x}+\frac{△x^2}{2}\frac{{\partial }^2(rU)}{\partial x^2}$$

$$+[1-\ell -r]U+\mathcal{O}\left(△x^3\right).$$

$$△t\frac{\partial U}{\partial t}+\hat{\beta }(x)$$

$$+\frac{△x^2}{2}\frac{{\partial }^2(\ell +r)U}{\partial x^2}+\mathcal{O}\left(△x^3,△t^2,\hat{\beta }(x{)}^2\right).$$

$${\color[rgb]{0,0,0}\frac{\partial U}{\partial t}+\beta(x)\frac{\partial U}{\partial\xi}\normalcolor}=\frac{\triangle x}{\triangle t}\frac{\partial(\ell-r)U}{\partial x}+\frac{\triangle x^{2}}{2\triangle t}\frac{\partial^{2}(\ell+r)U}{\partial x^{2}}+\mathcal{O}\left(\triangle x^{3},\triangle t,\hat{\beta}(x)^{2}\right).$$

$$U_t+\beta (x)U_{\xi }=\frac{△x^2}{2△t}{U}_{xx}+\mathcal{O}\left(△x^3,△t,\hat{\beta }(x{)}^2\right).$$

$$U_t+\beta (x)U_{\xi }=D{U}_{xx},$$

$$\{\begin{array}{ll}& U_t+\beta (x)U_{\xi }=D{U}_{xx},(x,\xi )\in {\Omega }^1,t>0\\ & U=\varphi (x,\xi ),(x,\xi )\in {\Omega }^1,t=0\\ & {\lim }_{|x|\to \infty }U=0,(x,\xi )\in {\Omega }^1,t>0.\end{array}$$

$$\{\begin{array}{ll}& U_t+\beta (\mathbf{x})U_{\xi }=D_n{\nabla }^{2}U,(\mathbf{x},\xi )\in {\Omega }^n,t>0\\ & U=\varphi (\mathbf{x},\xi ),(\mathbf{x},\xi )\in {\Omega }^n,t=0\\ & {\lim }_{|\mathbf{x}|\to \infty }U=0,(\mathbf{x},\xi )\in {\Omega }^n,t>0,\end{array}$$

$$\begin{array}{rl}V(x,t+& △t,\xi +\hat{\beta }(x))=\ell (x+△x)V(x+△x,t,\xi )\\ & +r(x-△x)V(x-△x,t,\xi )+[1-r(x)-\ell (x)]V(x,t,\xi )\\ & -H(\xi -{\xi }_c)V(x,t,\xi ),\end{array}$$

$$\{\begin{array}{ll}& V_t+\beta (x)V_{\xi }=D{V}_{xx}-\tilde{H}(\xi -{\xi }_c)V,(x,\xi )\in {\Omega }^1,t>0\\ & V=\varphi (x,\xi ),(x,\xi )\in {\Omega }^1,t=0\\ & {\lim }_{|x|\to \infty }V=0,(x,\xi )\in {\Omega }^1,t>0,\end{array}$$

$$p(x,t)={\int }_0^{\infty }V(\mathbf{x},t,\xi )d\xi .$$

$$p(\mathbf{x},t)={\int }_0^{{\xi }_c}U(\mathbf{x},t,\xi )d\xi .$$

$$\mathbb{P}(\mathbf{x},t)=\frac{{\int }_0^{{\xi }_c}U(\mathbf{x},t,\xi )d\xi }{{\int }_0^{\infty }U(\mathbf{x},t,\xi )d\xi }=\frac{p(\mathbf{x},t)}{{\int }_0^{\infty }U(\mathbf{x},t,\xi )d\xi }.$$

$${\int }_0^{{\xi }_c}{U}_{t}d\xi +{\int }_0^{{\xi }_c}\beta (\mathbf{x})U_{\xi }d\xi ={\int }_0^{{\xi }_c}{D}_n{\nabla }^{2}Ud\xi .$$

$$\{\begin{array}{ll}& p_t-D_n{\nabla }^{2}p=f(\mathbf{x},t),\mathbf{x}\in {\mathbb{R}}^n,t>0\\ & p(\mathbf{x},0)=g(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^n,t=0\\ & {\lim }_{|\mathbf{x}|\to \infty }p=0,\mathbf{x}\in {\mathbb{R}}^n,t>0,\end{array}$$

$$F(\mathbf{x},t)=-\frac{\partial }{\partial t}{\int }_{{\mathbb{R}}^n}p(\mathbf{x},t)d\mathbf{x}.$$

$$M={\int }_0^{\infty }tF(\mathbf{x},t)dt.$$

$$M={\int }_0^{\infty }{\int }_{{\mathbb{R}}^n}p(\mathbf{x},t)dxdt.$$

$$\frac{dE}{dt}$$

$$=-{\int }_{{\Omega }^n}\beta (\mathbf{x})U{U}_{\xi }d\xi d\mathbf{x}+D_n{\int }_{{\Omega }^n}U{\nabla }^{2}Ud\xi d\mathbf{x}.$$

$$\int_{\Omega^{n}}\beta(\mathbf{x})UU_{\xi}\,d\xi\,d\mathbf{x}=\int_{\mathbb{R}^{n}}\left[\beta(\mathbf{x})U^{2}\Big{|}_{\xi=0}^{\infty}\right]\,d\mathbf{x}-\int_{\Omega^{n}}\beta(\mathbf{x})UU_{\xi}\,d\xi\,d\mathbf{x}.$$

$${\int }_{{\mathbb{R}}^n}|\nabla U{|}^{2}d\mathbf{x}={\int }_{\partial {\mathbb{R}}^n}U\nabla U\cdot \hat{\eta }d\mathbf{x}-{\int }_{{\mathbb{R}}^n}U{\nabla }^{2}Ud\mathbf{x},$$

$${\int }_{{\Omega }^n}U{\nabla }^{2}Ud\xi d\mathbf{x}=-{\int }_{{\Omega }^n}|\nabla U{|}^{2}d\xi d\mathbf{x}.$$

$$\frac{dE}{dt}=-D_n{\int }_{{\Omega }^n}|\nabla U{|}^{2}d\xi d\mathbf{x}\le 0.$$

$$\frac{\partial }{\partial t}{\int }_{{\Omega }^n}Ud\xi d\mathbf{x}$$

Full text

KEYWORDS: Agent-based models, Upscaling, Green’s functions, Operator splitting, Chemical absorption

1 Introduction

In many applications, we may wish to track the individual movement of a collection of cells or agents that exhibit state changes and interact with a dynamic environment. A Lagrangian framework such as a cellular automaton or an agent-based model can be used to track the movement of each cell, while also accounting for state changes [11, 12, 14, 20]. To implement this type of model, one would define the cell as the agent and a new location at each time point is chosen based on pre-defined rules. With current computer architectures, simulating millions of cells is feasible but is offset by the computational time required to compute a sufficient number of simulations for analysis [1, 8, 24, 29]. In addition, in the case where one wants to couple movement rules to evolving profiles on the domain, as well as tracking cumulative variables related to environmental conditions in multiple dimensions, computational time could be prohibitive.

Several other modeling approaches can also be utilized to understand the dynamics of a large group of cells or agents that exhibit state changes and interact with their environment; each of these methods has their own limitations and advantages. The transition between cell states can be modeled as a stochastic process, and in the case of a Markovian or memoryless process, the transition rate will only depend on the current state, following a Poisson distributed random process [13, 23]. In terms of a discrete-time Markov Chain, each state can correspond to a particular combination of cell state and cell location. Using this type of approach, the probability of a cell being in a given state at a given time can be determined and the master equation for the rate of change of the associated probability can be obtained. Since analysis of cellular processes is often times easier in the continuous setting, one can easily move from discrete probabilities of cell states to a continuous probability density of cell states by assuming a continuum of cell states. However, in order to keep track of a cumulative variable accounting for interactions with the environment, this approach would break down if memory in the system was necessary.

Random walk models are also stochastic processes that consist of sequential random steps of movement; they have been widely used to investigate cellular motility, often in a spatially homogeneous environment [3, 7, 15, 30, 31]. Assuming the moving agent or cell is memoryless, an equation governing the spatiotemporal evolution of the density of cells can be determined. It corresponds to a standard diffusion equation if there is no bias in the motion or an advection-diffusion equation if there is bias in the motion [3, 7, 15, 31]. The advection-diffusion equation can capture different taxis, biasing the probability of movement based on chemical profiles (chemotaxis), temperature gradients (thermotaxis), fluid flow (rheotaxis), or environmental mechanical stiffness (durotaxis) [2, 6, 17, 18, 19]. The continuum limit of the stochastic process is often formulated in the case of cell motility since it is tractable from an analytical perspective and we have existing computational methods to easily solve these governing equations. In this framework, accounting for different cell states would correspond to a system of coupled partial differential equations where local sinks or sources would describe leaving one cell state and entering another cell state. Currently, there is no random walk modeling framework to account for cells changing states due to a cumulative environmental coupling.

Cellular processes, such as absorbing chemicals or nutrients, moving, and state transitioning will depend on the local and dynamic environment [16, 22, 25, 26]. Often, exposure to a chemical or drug beyond a threshold will cause a change in state or motion, thus it is important to capture the cumulative chemical exposure. To date, analysis has primarily focused on the motion of a single cell type in a homogeneous environment or a sequence of cell states in time (not accounting for motion) [7, 23]. On the other hand, stochastic agent-based models can simulate many cells and states, but analyzing these models is intractable and computations can become prohibitive with a large number of cells coupled to the evolving environment. Hence, one approach to resolve this issue is by implementing multiscale models that bridge the gap between the individual cell level and dynamics at the global continuum level. For example, recent work has focused on developing movement rules coupled to chemical concentrations in the case of bacterial chemotaxis [5, 34, 33].

In this paper, we focus on the development of new computational methods to capture the continuum density of cells or agents in space, accounting for the cumulative exposure to the environment as a continuous variable. As outlined in Section 2, our motivating example will be a cell that moves and absorbs chemicals; a state change occurs when the cell has absorbed a critical (toxic) threshold of chemicals, causing the cell to die. In Section 3, we will show how the new governing equation is able to capture these dynamics. An analysis of this equation is detailed in Section 4 and the numerical method is outlined in Section 5. Representative numerical results are given in Section 6, comparing computation of our new governing equations with the corresponding agent-based model for the case of cells that randomly move and absorb chemical from the surrounding environment. Additional commentary is given for several limiting cases in Section 7.

Notation. $\mathbb{N}$ denotes the natural numbers and $\mathbb{R}^{n}$ denotes $n$-dimensional Euclidean space. All vectors will be denoted with bold face, i.e. $\mathbf{x}=[x_{1},x_{2},\ldots,x_{n}]^{T}$ when $\mathbf{x}\in\mathbb{R}^{n}$ and the superscript $T$ denotes the vector transpose. The $L^{p}$ space is defined by generalizing the $p$-norm for vector spaces $\mathbb{R}^{n}$, whereas $C^{k}$ defines the space of continuously differentiable functions up to the $k^{th}$ derivative. For ease of notation, we define $\Omega^{n}=\mathbb{R}^{n}\times[0,\infty)$, the spatial and absorption domains. The convolution of two functions $f$ and $g$ is denoted as $f*g$ and we define $*^{m}$ as

[TABLE]

We will use $\delta(\cdot)$ to denote the Dirac delta distribution. To simplify, we will denote $||\cdot||_{1}\equiv\int_{\Omega^{n}}|\cdot|\,d\xi\,d\mathbf{x}$ and $||\cdot||_{2}\equiv\int_{\Omega^{n}}(\cdot)^{2}\,d\xi\,d\mathbf{x}$. The error function is defined as $erf(z)\equiv(2/\sqrt{\pi})\int_{0}^{z}e^{-t^{2}}dt$. The evaluation of $\lfloor d\rfloor$ gives the greatest integer that is less than or equal to $d$.

2 Problem formulation and preliminaries

Suppose a spatial domain contains a spatially-varying (or time-varying) chemical concentration. For simplicity, we will assume the chemical concentration is a positive, spatially-dependent distribution, $C(\mathbf{x})$. We then insert a cell in this domain at location $\mathbf{x}_{0}$ - to avoid semantic confusion later in this paper, we will refer to this cell as an agent. A schematic of this setup is shown in Fig. 1. This agent has a given probability of moving at each time point. In the 2-dimensional setup, the agent may remain stationary or move left, right, up, or down as shown with the dotted arrows in Fig. 1. After moving, the agent will absorb a certain amount of chemical according to a function $\hat{\beta}(x)$, which depends on the local chemical concentration. Although we are not fixing a specific form of the function $\hat{\beta}(x)$ for our analysis, we will assert that this function preserves the property that if $C(x)>0$, then $\hat{\beta}(x)>0$. When the cumulative absorption within the agent reaches a critical threshold, the agent changes state (e.g. the agent dies).

At first glance, we may consider the agent in our absorption model as having only two states: live and dead. However, it is subtly more complicated. Considering that the chemical concentration could be spatially and/or time-dependent, the particular path along which an agent travels may affect the amount of chemical the agent absorbs. For example, suppose there is no chemical concentration to the left of $\mathbf{x}_{0}$ and there is chemical to the right of $\mathbf{x}_{0}$. Further, consider two distinct paths the agent may travel: in one path, the agent is contained within the left side of the domain and terminates at $\mathbf{x}_{0}$ at time $t$, whereas in another path, the agent is contained within the right side of the domain and terminates at $\mathbf{x}_{0}$ at $t$. The agent which traveled along the first path will not have absorbed any chemical by time $t$, but the agent which traveled along the second path will absorb some chemical particles.

We need to account for varying amounts of chemical concentration by treating the amount absorbed as distinct states. However, as opposed to a compartmental model, the chemical concentration is a continuous variable. In order to account for this, we need to consider cumulative absorption as a dimension orthogonal to both the temporal and spatial dimensions. That is, the agent is at a location $\mathbf{x}$, having a cumulative chemical absorption $\xi$, at a particular time $t$.

The cumulative amount absorbed is path dependent. However, we cannot say that $\xi$ is dependent on space or time, just as we cannot say that $\mathbf{x}$ is dependent on time. All three variables are linked by our model, but should be considered independent.

3 The continuum model

Through upscaling, we will derive the continuum absorption model. We develop a discrete difference equation from the agent-based model (ABM), and then take the continuum limit, an approach used for standard random walk models [7]. At each iteration, our ABM agent $\mathfrak{u}$ moves in the domain with spatial-step $\triangle x$ and time-step $\triangle t$. Then, $\mathfrak{u}$ absorbs chemical particles at its new location $\mathbf{x}$ based on a spatially-dependent function $\hat{\beta}(\mathbf{x})$ that depends on the chemical distribution, $C(\mathbf{x})$. We interpret $\hat{\beta}(\mathbf{x})$ as the amount the cell absorbs during the time step of length $\triangle t$.

3.1 Derivation of single-state absorption model

Assume the ABM agent $\mathfrak{u}$ is initialized in the 1-D spatial domain $\mathbb{R}$. We define the variable $U(x,t,\xi)$ to denote the density of agent $\mathfrak{u}$ at location $x$ at time $t$, having absorbed $\xi$ total particles. Suppose $\mathfrak{u}$ moves a distance $\triangle x$ to the right to be at $x$ at the next time increment, $t+\triangle t$. We then have that $\mathfrak{u}$ absorbs $\hat{\beta}(x)$ chemical particles at $x$, with a cumulative absorption of $\xi+\hat{\beta}(x)$ particles at time $t+\triangle t$ (when $\mathfrak{u}$ had a cumulative absorption of $\xi$ at time $t$). Similarly, if $\mathfrak{u}$ moves a distance $\triangle x$ to the left to be at $x$ at the next time increment, then $\mathfrak{u}$ will have cumulatively absorbed $\xi+\hat{\beta}(x)$ particles. Letting $\ell(x)$ be the probability of moving left at location $x$ and $r(x)$ the probability of moving right at $x$, we can assert a difference equation modeling this behavior:

[TABLE]

where $\ell(x)+r(x)\leq 1$ for all $x\in\mathbb{R}$. The $U(x,t+\triangle t,\xi+\hat{\beta}(x))$ term expresses the fact that agent $\mathfrak{u}$ absorbs an additional $\hat{\beta}(x)$ chemical particles at location $x$ and time $t+\triangle t$. The right-hand side of (1) accounts for all the different possible ways (along with their respective probabilities) that agent $\mathfrak{u}$, having absorbed $\xi$ total number of chemical particles at time $t$ can be at location $x$ at time $t+\triangle t$.

Assuming $U\in C^{2}(\mathbb{R},[0,\infty),[0,\infty))$ and $q(x)\in C^{2}([0,1])$ for $q(x)$ as either $\ell(x)$ or $r(x)$, we can perform a Taylor expansion on $U$ and the moving probability $q(x)$ in (1) and get

[TABLE]

Inserting these expansions into (1) results in

[TABLE]

Rearranging terms and simplifying gives us

[TABLE]

If we define $\beta(x)=\hat{\beta}(x)/\triangle t$ and rearrange terms, we have

[TABLE]

For simplification, we will let $\ell(x)=r(x)=1/2$ for all $x\in\mathbb{R}$. Our equation then reduces to:

[TABLE]

We assume that $\triangle t\sim\triangle x^{2}$ and $\hat{\beta}(x)=\mathcal{O}(\triangle t)$. Taking the limit as $\triangle t,\triangle x\to 0$ results in the following governing continuum equation:

[TABLE]

where $D=\lim_{\triangle t,\triangle x\to 0}\triangle x^{2}/(2\triangle t)$. For this paper, we assume far-field boundary conditions and the initial condition can depend on $x$ and $\xi$.

Hence, our partial differential equation (PDE) for chemical absorption is as follows

[TABLE]

In a similar way, assuming that the spatial step $\triangle x$ is the same in every direction, we can derive a continuum PDE in $n$ spatial dimensions. The resulting PDE is as follows:

[TABLE]

where $D_{n}=\lim_{\triangle x,\triangle t\to 0}\triangle x^{2}/(2n\triangle t)$. Here we should note that the above equations (5) and (6) only model cumulative absorption, without taking into account any state transitions due to this cumulative absorption. We address two possible methods to account for state changes in the following subsection.

3.2 Absorption threshold

Now that we have a governing equation for tracking the cumulative absorption property of an agent, we can address possible absorption-dependent state-changes. Suppose the agent changes state if the cumulative chemical absorption is greater than some absorption capacitance, $\xi_{c}$. That is, a cell is initially in the live state if $\xi<\xi_{c}$ and switches to a different state, possibly dying, if $\xi\geq\xi_{c}$. We will denote the spatially-varying, total density of agents in the live state as $p(\mathbf{x},t)$.

We have two possible methods111We consider methods that keep $\beta$ from being dependent on the cumulative absorption amount, $\xi$. One could derive a formula where $\beta$ depends on both the chemical concentration, $C(x)$, and $\xi$ to ensure that cells do not absorb chemical beyond the threshold. However, the derivation and resulting Taylor series would produce additional terms that could make the analysis and numerical solutions more difficult.

for calculating the density of cells in a live state: by modifying our PDE model or by solving the PDE model and integrating $\xi$ over the domain of interest, $[0,\xi_{c}]$. First, we will examine modifying our PDE model and in this case, we will use the variable $V(x,t,\xi)$ to denote the density of agents in the live state (rather than the variable $U$ to avoid confusion in later sections). Our difference equation includes a Heaviside function centered at $\xi_{c}$, since the agent switches state if $\xi\geq\xi_{c}$,

[TABLE]

where $H(\xi-\xi_{c})=\begin{cases}1&:\text{ if }\xi\geq\xi_{c}\\ 0&:\text{ otherwise}\end{cases}$ denotes the Heaviside function centered at $\xi_{c}$. With the same assumptions as in the previous section and expanding in a Taylor series, we obtain the following PDE modeling live agents following an unbiased random walk:

[TABLE]

where $\tilde{H}(\xi-\xi_{c})=\lim_{\triangle t\to 0}\frac{1}{\triangle t}H(\xi-\xi_{c})=\begin{cases}+\infty&:\text{ if }\xi\geq\xi_{c}\\ 0&:\text{ otherwise}\end{cases}$ is the result of taking the continuum limit. We can then calculate the spatially-variable, total density of cells in the live state as

[TABLE]

This method offers us the ability to develop more complex models, such as when $\beta$ can be negative, when the transition is stochastic (in which case we would substitute the Heaviside function with a probability distribution), or when we are interested in modeling agents in the secondary state. In this case, we can more easily separate different populations and capture additional state change driven dynamics (e.g. where movement rules could depend on the given state).

The alternate approach, which is the focus of this paper, is suitable in the case where we are only interested in the live agents and there is a single deterministic state transition (when $\xi>\xi_{c}$) and $\beta>0$. In this case, it is possible to calculate the spatially-variable, total density, $p$, directly from the absorption model, $U$ as

[TABLE]

If we initialize $\int_{\Omega^{n}}U(\mathbf{x},0,\xi)\,d\xi\,d\mathbf{x}=1$, then we can consider $p(\mathbf{x},t)$ the probability that an agent is at location $\mathbf{x}$ at time $t$ and in the initial live state222If we want to find the probability an agent is at a particular location at a given time, given that the agent is in the initial live state, we can calculate:

$\mathbb{P}(\mathbf{x},t)=\frac{\int_{0}^{\xi_{c}}U(\mathbf{x},t,\xi)\,d\xi}{\int_{0}^{\infty}U(\mathbf{x},t,\xi)\,d\xi}=\frac{p(\mathbf{x},t)}{\int_{0}^{\infty}U(\mathbf{x},t,\xi)\,d\xi}.$

.

We can rewrite (6) as a PDE of $p$. Let us integrate the terms from 0 to $\xi_{c}$ with respect to $\xi$. This gives us

[TABLE]

Given $U\in L^{1}(\Omega^{n})$, we can switch derivatives and integrals using Fubini’s theorem. The above system reduces to a non-homogeneous diffusion equation

[TABLE]

where f(\mathbf{x},t)=-\beta(\mathbf{x})U\Big{|}_{\xi=\xi_{c}} and $g(\mathbf{x})=\int_{0}^{\xi_{c}}\phi(\mathbf{x},\xi)$. If we know the value of $f$, then we have an explicit solution for $p$ using the method of Green’s functions (fundamental solutions). In most cases we will not have the explicit value of $f(\mathbf{x},t)$, in which case we must first solve for $U(\mathbf{x},t,\xi)$ before integrating to compute $p(\mathbf{x},t)$.

We can use the value of $p$ to calculate cellular properties of interest, such as flux out of the initial live state or the average time in the initial live state.

3.3 Mean occupancy time

We may be interested in the mean time an agent is in the initial live state, which is denoted as the mean occupancy time (MOT). In a manner similar to deriving the mean first passage time, this is the first moment of the total flux out of a particular state.

The total flux out of the initial state can be computed as

[TABLE]

The negative sign is due to the fact that we are tracking the density exiting the initial state. It follows that the MOT is

[TABLE]

Since $p\in L^{1}(\mathbb{R}^{n})$ and for any finite location $\mathbf{x}\in\mathbb{R}^{n}$, $\lim_{t\to\infty}p(\mathbf{x},t)=0$, we can use integration by parts to derive the MOT,

[TABLE]

4 Mathematical analysis

Through the derivation of this continuous approximation, higher order terms in the Taylor series expansions were neglected. We must still ensure that we are maintaining the proper physics with this new equation. For example, we wish that energy in the system is not increasing and that the total quantity of agents or cells is conserved. In addition, since the governing equation (6) is classified as a mixed Parabolic-Hyperbolic PDE, there is no generalized theorem we can apply to show it is well-posed. To this end, Theorems 3-5 in this section will prove existence, uniqueness, and continuous dependence on initial data, respectively.

4.1 Energy & conservation

In order to prove uniqueness and the continuous dependence of the PDE solution on initial data, we need to show that there is some time-dependent functional $E(t)$, such that our solution $U$ of (6) satisfies $0\leq E(t)\leq E(0)$ for all $t>0$. We will refer to this functional as the energy of the solution at time $t$. To match the physics of the ABM simulation, the energy of our PDE should be non-increasing, but we need to prove that our PDE does not lose this feature during the process of deriving the continuum approximation.

Theorem 1**.**

Suppose $\beta(\mathbf{x})>0$ for all $\mathbf{x}\in\mathbb{R}^{n}$. The PDE (6) with the energy functional 333The energy functional is not meant to be interpreted as physical energy (i.e. it is not kinetic or potential energy), it is a naming convention that is consistent with classical PDE theory [9]. $E(t)=\frac{1}{2}||U||_{2}$ satisfies the inequality $0\leq E(t)\leq E(0)$.

Proof.

Via a calculation,

[TABLE]

First, integration by parts in the variable $\xi$ gives

[TABLE]

We assume that for any finite $t>0$ that $U=0$ as $\xi\to\infty$. Given $\beta(\mathbf{x})>0$, then $U=0$ at $\xi=0$ for any $t>0$. Thus, $\int_{\Omega^{n}}\beta(\mathbf{x})UU_{\xi}\,d\xi\,d\mathbf{x}=0$. Second, by the Divergence product rule,

[TABLE]

where $\hat{\eta}$ is the unit outward normal vector. Considering $\lim_{|\mathbf{x}|\to\infty}U=0$, we have that

[TABLE]

Therefore, we have that for every $t>0$,

[TABLE]

Seeing that $\frac{dE}{dt}\leq 0$, we have $0\leq E(t)\leq E(0)$.

∎

Note that the energy functional $E_{V}(t)=\frac{1}{2}||V||_{2}$ also satisfies the inequality $0\leq E_{V}(t)\leq E_{V}(0)$ since for every $t>0$, $\int_{\Omega^{n}}\tilde{H}(\xi-\xi_{c})U^{2}\,d\xi\,d\mathbf{x}\geq 0$.

In the ABM simulation, no agent is removed from the system. Again, we want the PDE solution to match the important physics of the ABM simulation. We do so by proving that the solution $U$ is conserved at each time $t$ over the entire domain $\Omega^{n}$.

Theorem 2**.**

(Conservation) Suppose $U\in L^{1}(\Omega^{n})$ solves (6) and $\beta(\mathbf{x})>0$ for all $\mathbf{x}\in\mathbb{R}^{n}$. Then $\int_{\Omega^{n}}U\,d\xi\,d\mathbf{x}=\int_{\Omega^{n}}\phi(\mathbf{x},\xi)\,d\xi\,d\mathbf{x}$ for any $t>0$.

Proof.

By means of a calculation,

[TABLE]

Since $\lim_{|\mathbf{x}|\to\infty}U=0$ we have that the first term is 0. Also, given $\beta(\mathbf{x})>0$ for all $x\in\mathbb{R}^{n}$, then for $t>0$ we have $U(\mathbf{x},t,\xi=0)=0$, and the second term is also 0. It follows that $\frac{\partial}{\partial t}\int_{\Omega^{n}}U\,d\xi\,d\mathbf{x}=0$. Therefore, $U$ is conserved. ∎

4.2 Operator-splitting semi-discrete solution

We will approximate a solution to the PDE in (6) by splitting the linear operator and then solving the resulting system iteratively. This gives us a solution that is discrete in time and continuous in spatial and absorption dimensions. We will first derive this semi-discrete solution and then show that it is well-posed.

Let $U=\widehat{U}(\mathbf{x},t|\xi)\overline{U}(\xi,t|\mathbf{x})$, where $\widehat{U}$ leaves $\xi$ fixed and $\overline{U}$ leaves $\mathbf{x}$ fixed. We can see that $\overline{U}\widehat{U}_{t}+\widehat{U}\overline{U}_{t}+\beta(\mathbf{x})\widehat{U}\overline{U}_{\xi}=D_{n}\overline{U}\nabla^{2}\widehat{U}$ and it follows that $\overline{U}\left(\widehat{U}_{t}-D_{n}\nabla^{2}\widehat{U}\right)+\widehat{U}\left(\overline{U}_{t}+\beta(\mathbf{x})\overline{U}_{\xi}\right)=0.$ Assuming that $\overline{U}$ and $\widehat{U}$ are not identically 0, we can then solve the following PDEs

[TABLE]

[TABLE]

We solve the system in (14) using the method of Green’s functions444 The Green’s function used in this paper is the fundamental solution of the diffusion equation. However, for other spatial domains, such as a half plane or disk, we can use the method of images with the fundamental solution to derive an appropriate Green’s function. and convoluting with the initial condition:

[TABLE]

where

[TABLE]

is the fundamental solution of the diffusion equation in $\mathbb{R}^{n}$. We solve (15) using the method of characteristics:

[TABLE]

Our solution of (6) alternates between (16) and (18) as the solution marches forward in time. As we are not solving the system simultaneously, we will choose a length of time, $0<\tau\ll 1$, in which each solution is valid. We will denote the solution at time $t=m\tau$ as $U^{m}(\mathbf{x},\xi)$. The following iterative algorithm solves the semi-discrete, operator splitting system:

• •

Initialize $U^{0}(\mathbf{x},\xi)=\phi(\mathbf{x},\xi)$

• •

For $m=1,2,\ldots$:

• $\lozenge$

$\overline{U}^{m-1}(\mathbf{x}|\xi)=U^{m-1}(\mathbf{x},\xi)$

• $\lozenge$

$\widehat{U}^{m}(\xi|\mathbf{x})=\overline{U}^{m-1}(\mathbf{x}|\xi-\beta(\mathbf{x})\tau)$

• $\lozenge$

$U^{m}(\mathbf{x},\xi)=G(\mathbf{x},\tau)*\widehat{U}^{m}(\xi|\mathbf{x})$

Combining these solutions gives us the recurrence relation for the semi-discrete solution with time step $\tau$:

[TABLE]

Additionally, we can use our recurrence relation to rewrite the solution at $t=m\tau$ in terms of the initial condition $\phi(\mathbf{x},\xi)$, given as

[TABLE]

4.3 Existence

Our semi-discrete solution for $U^{m}(\mathbf{x},\xi)$, given in (20), depends on the recurrence time step $\tau$ and the number of iterations $m$. So, we define

[TABLE]

as the approximation555 We could also make a similar proof for approximating V with $z_{m,\tau}(\mathbf{x},\xi)\equiv G(\mathbf{x},\tau)*^{m}\left[\phi(\mathbf{x},\xi-\beta(\mathbf{x})m\tau)\mathbf{1}_{[0,\xi_{c})}(\xi-\beta(\mathbf{x})m\tau-\xi_{c})\right]$. The indicator function is a result of solving the same operator splitting system as U with the additional equation

$\begin{cases}\tilde{V}_{t}=-\tilde{H}(\xi-\xi_{c})\tilde{V},\quad&\mathbf{x}\in\mathbb{R}^{n},\xi\in[0,\infty),t>0\\ \tilde{V}=\bar{V},\quad&\mathbf{x}\in\mathbb{R}^{n},\xi\in[0,\infty),t=0.\end{cases}$

of $U(\mathbf{x},m\tau,\xi)$, accounting for the choices of both $\tau$ and $m$ (the same notation as used in [4]). We want to show that the $L^{1}(\Omega^{n})$ limit, $\mathcal{U}(\mathbf{x},t,\xi)$, of the sequence $\left\{z_{m,t/m}\right\}_{m\in\mathbb{N}}$ exists. That is, the limit of the recurrence relation for a time $t$ exists when the recurrence time step, $\tau=t/m$, approaches 0. If this limit exists, it proves the existence of a solution, $U(\mathbf{x},t,\xi)$, to the governing PDE, given in (6).

For the $L^{1}$ limit to make sense, we first need to show that $z_{m,\tau}\in L^{1}(\Omega^{n})$.

Lemma 1**.**

Suppose $\phi\in L^{1}(\Omega^{n})$ and $z_{m,\tau}(\mathbf{x},\xi)=G(\mathbf{x},\tau)*^{m}\phi(\mathbf{x},\xi-\beta(\mathbf{x})m\tau)$ for all $m\in\mathbb{N},\tau>0$. Then $z_{m,\tau}(\mathbf{x},\xi)\in L^{1}(\Omega^{n})$, for all $\mathbf{x}\in\mathbb{R}^{n},\xi\geq 0$.

Proof.

We know that $G(\mathbf{x},\tau)\in L^{1}(\Omega^{n})$. By reason that $L^{1}(\Omega^{n})$ is closed under convolution, we have that $G(\mathbf{x},\tau)*^{m}\phi(\mathbf{x},\xi-\beta(\mathbf{x})m\tau)\in L^{1}(\Omega^{n})$.

∎

Lemma 2**.**

For any $t>0$, $G(\mathbf{x},t/m)*^{m}\delta(\xi)\to\delta(\mathbf{x},\xi)$ as $m\to\infty$ in $L^{1}(\Omega^{n})$.

Proof.

We know that $\lim_{m\to\infty}\int_{\mathbb{R}^{n}}\left|G(\mathbf{x},t/m)-\delta(\mathbf{x})\right|\,d\mathbf{x}=0$, so it follows that $\lim_{m\to\infty}||G(\mathbf{x},t/m)*\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$. First, we will show, via induction, that $\lim_{m\to\infty}||G(\mathbf{x},t/m)*^{h}\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$ for any $h\in\mathbb{N}$.

As a base case, we will show that $\lim_{m\to\infty}||(G*G)(\mathbf{x},t/m)*\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$. By the Dominated Convergence theorem, we can see that

[TABLE]

By calculation,

[TABLE]

Now, we can assume that $\lim_{m\to\infty}||G(\mathbf{x},t/m)*^{h}\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$ for some $h\in\mathbb{N}$. We can calculate that

[TABLE]

by the inductive assumption. It follows that

[TABLE]

for all $h\in\mathbb{N}$.

If we define $h=m$ for $m\in\mathbb{N}$. Then $\lim_{m\to\infty}||G(\mathbf{x},t/m)*^{m}\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$. ∎

Lemma 3**.**

Suppose $\phi\in L^{1}(\Omega^{n})$ and $z_{m,t/m}(\mathbf{x},\xi)=G(\mathbf{x},t/m)*^{m}\phi(\mathbf{x},\xi-\beta(\mathbf{x})t)$ for any $m\in\mathbb{N},t>0$, as defined in (21). Then $\left\{z_{m,t/m}\right\}_{m\in\mathbb{N}}$ is a Cauchy sequence in $L^{1}(\Omega^{n})$.

Proof.

We want to show

[TABLE]

From Lemma 2 we know that, for any $m\in\mathbb{N}$, $\lim_{m\to\infty}||G(\mathbf{x},t/m)*^{m}\delta(\xi)-\delta(\mathbf{x},\xi)||_{1}=0$. Then

[TABLE]

It follows that

[TABLE]

Therefore, $\lim_{p,q\to\infty}||z_{p,t/p}(\mathbf{x},\xi)-z_{q,t/q}(\mathbf{x},\xi)||_{1}=0$. ∎

Theorem 3**.**

(Existence) Suppose $\phi\in L^{1}(\Omega^{n})$. There exists a solution, $U\in L^{1}(\Omega^{n})$, to the governing PDE:

[TABLE]

Proof.

Choose any $(\mathbf{x},\xi)\in\Omega^{n}$ and any $t>0$. Suppose $\phi\in L^{1}(\Omega^{n})$ and define $z_{m,\tau}(\mathbf{x},\xi)=G(\mathbf{x},\tau)*^{m}\phi(\mathbf{x},\xi-\beta(\mathbf{x})m\tau)$ for any $m\in\mathbb{N}$, $\tau>0$. By Lemma 1 and Lemma 3 we know that $z_{m,t/m}\in L^{1}(\Omega^{n})$ is a Cauchy sequence. On account of $L^{1}$ being complete, there exists a $\mathcal{U}\in L^{1}(\Omega^{n})$ such that $\lim_{m\to\infty}||z_{m,t/m}(\mathbf{x},\xi)-\mathcal{U}(\mathbf{x},t,\xi)||_{1}=0$. Since $z_{m,t/m}(\mathbf{x},\xi)$ satisfies the operator-split PDE for all $m\in\mathbb{N}$, we know that $\mathcal{U}(\mathbf{x},t,\xi)$ satisfies the operator-split PDE. Therefore $\mathcal{U}(\mathbf{x},t,\xi)$ satisfies the time-continuous PDE. ∎

4.4 Uniqueness & continuous dependence on initial data

Theorem 4**.**

(Uniqueness) The solution to PDE (6) is unique.

Proof.

Suppose we have two solutions, $U_{1},U_{2}\in L^{1}(\Omega^{n})$ to the PDE (6). We will define $W=U_{2}-U_{1}$. Given (6) is linear, we know $W(\mathbf{x},t,\xi)$ solves the PDE:

[TABLE]

From the energy argument in Theorem 3, we know that

[TABLE]

Seeing that $E_{w}(0)=\frac{1}{2}\int_{\Omega^{n}}W(\mathbf{x},0,\xi)^{2}\,d\xi\,d\mathbf{x}=0$, we know that $E_{w}(t)=0$ for all $t$. By definition of $E_{w}(t)$, we demonstrated that

[TABLE]

Therefore, $U_{1}(\mathbf{x},t,\xi)=U_{2}(\mathbf{x},t,\xi)$ almost everywhere. ∎

Theorem 5**.**

(Continuous Dependence on Initial Data) Consider any $\epsilon>0$. Suppose $U_{1}$ satisfies the PDE

[TABLE]

and $U_{2}$ satisfies the PDE

[TABLE]

where $||\phi_{1}(\mathbf{x},\xi)-\phi_{2}(\mathbf{x},\xi)||_{2}<\epsilon$. Then $||U_{1}-U_{2}||_{2}<\epsilon$.

Proof.

We will define $W=U_{1}-U_{2}$. As (24) and (25) are both linear, $W$ solves the PDE

[TABLE]

Let us define the energy of (26) as

[TABLE]

By the same argument in the proof of Theorem 3, we have that $0\leq||W(\mathbf{x},t,\xi)||_{2}\leq||W(\mathbf{x},0,\xi)||_{2}$. Since $||W(\mathbf{x},0,\xi)||_{2}=||\phi_{1}(\mathbf{x},\xi)-\phi_{2}(\mathbf{x},\xi)||_{2}<\epsilon$ and $||W(\mathbf{x},t,\xi)||_{2}=||U_{1}(\mathbf{x},,t,\xi)-U_{2}(\mathbf{x},t,\xi)||_{2}$, we have that

[TABLE]

∎

From the energy argument in Theorem 3, we know that the PDE for V also satisfies the properties of uniqueness and continuous dependence on initial data. Thus, the model for V is also well-posed.

5 Numerical approximation

5.1 Fully discrete derivation

We are primarily interested in calculating the spatially-variable, total density of agents in the live state, $p$. While demonstrating that the PDE models for $U$ and $V$ are both well-posed, we explained how the only difference in approximate solutions was multiplying by an indicator function $\mathbf{1}_{[0,\xi_{c})}(\xi)$. This leads to the fact that the value of $p(\mathbf{x},t)$ using either method will be identical. So for our numerical computation, we will compute the spatially-variable, total density of agents in the initial live state with $p(\mathbf{x},t)=\int_{0}^{\xi_{c}}U(\mathbf{x},t,\xi)\,d\xi$. We will derive this numerical approximation within the spatial domain in 1-D, but the method can easily extend to higher dimensions. Considering that we will first solve (6) for $U$, we will discretize the region $\Omega^{1}$ using cell volumes (as opposed to discrete nodes). We divide the spatial component into $N$ bins666 The analytic solution requires the spatial domain to be $\mathbb{R}$. However, numerically, we need to choose a finite domain. We choose $N$ such that $G(N\,dx/2,\tau)$ is bounded close to 0. Similarly, the absorption domain is $[0,\infty)$. Our solution of interest is within the domain $[0,\xi_{c})$ and hence, we choose $K$ such that $K\,d\xi$ is larger than $\xi_{c}$. of width $dx$ and the absorption component into $K$ bins of width $d\xi$; the cell volumes have area $dxd\xi$.

These cell volumes will be defined as $\omega_{i,k}=B(x_{i},dx/2)\times[\xi_{k},\xi_{k+1}]$, where $dx$ is the spatial discretization step-size and $B(x_{i},dx/2)=\{y\in\mathbb{R}:|x_{i}-y|<dx/2\}$. For the following derivations, the spatial location will be indexed by $i$ and the cumulative absorption amount will be indexed by $k$. We then define $u_{i,k}^{m}\approx U^{m}(x_{i},\xi_{k})$ as

[TABLE]

the average value of $U^{m}$ in the cell volume $\omega_{i,k}$. Note that the continuous and semi-discrete solution is capitalized, $U(x,t,\xi)$ or $U^{m}(x,\xi)$, whereas the fully discrete solution is in lower-case, $u_{i,k}^{m}$.

We know the semi-discrete recurrence relation $U^{m+1}(x,\xi)=G(x,\tau)*U^{m}(x,\xi-\beta(x)\tau)$ from equation (19). This solution is fully discretized by integrating over the cell volume $\omega_{i,k}$. By recalling that $u_{i,k}^{m}$ is piece-wise continuous over $\omega_{i,k}$, we can solve the convolution exactly with the approximated solution:

[TABLE]

Since $u_{i,k}^{m+1}=\frac{1}{dx\,d\xi}\int_{\omega_{i,k}}U^{m}(y,z)\,dy\,dz$, we have

[TABLE]

where $d_{i}=\lfloor k+\beta(x_{i})\tau\rfloor$, the new absorption index. By calculation, we find that

[TABLE]

Our numerical method is then

[TABLE]

We discretize the density $p(x,t)$ as

[TABLE]

where $\mathcal{A}=\{k:kd\xi<\xi_{c}\}$. Therefore, we can represent $p$ numerically as $p_{i}^{m}=\sum_{k\in\mathcal{A}}u_{i,k}^{m}$, the exact integral using our piece-wise constant approximate solutions.

5.2 Stability

We can break down the numerical method into two steps: a diffusive step where we perform the convolution,

[TABLE]

and an absorption step $u_{i,d_{i}}^{m+1}=v_{i,k}^{m+1}$ where we change the indexing. Further, we define the discrete energy functional as

[TABLE]

To prove stability, we want to show $\mathbf{E}^{m+1}-\mathbf{E}^{m}\leq 0$. Note that due to the indexing change, for any $i$, \sum_{k=0}^{K}\left(u_{i,k}^{m+1}\right)^{2}\leq\sum_{k{\color[rgb]{0,0,0}=\normalcolor}0}^{K}\left(v_{i,k}^{m+1}\right)^{2}.

We can rewrite our numerical scheme as a matrix-vector product, $\mathbf{v}_{k}^{m+1}=\mathbf{u}_{k}^{m}*G=A\mathbf{u}_{k}^{m}$, where our discrete convolution matrix $A$ and vector indexing of $\mathbf{u}_{k}^{m}$ are the following

[TABLE]

given our definition of $G_{j}$, as defined in (28).

The difference between the energy functional at subsequent times, having absorbed $\xi\in[kd\xi,(k+1)d\xi]$ particles, is:

[TABLE]

where $E_{k}^{m}\equiv\sum_{i=0}^{N-1}\left(u_{i,k}^{m}\right)^{2}$.

Theorem 6**.**

The spectrum of $A^{T}A-I$ is {\color[rgb]{0,0,0}s\normalcolor}\left(A^{T}A-I\right)\leq 0, with the matrix $A$ defined in (31) and the scalars $G_{j}$ defined in (28).

Proof.

The discrete convolution matrix $A$ is a circulant matrix, so it has eigenvalues

[TABLE]

for $j=0,1,\ldots,N_{1}$, where $\gamma_{j}=\exp\left\{\frac{2\pi j}{N}\sqrt{-1}\right\}$ (the $N$-th root of unity). It follows that the amplitude of the $j$-th eigenvalue is

[TABLE]

Given $G_{\ell}\geq 0$ for all $\ell$, we have that $|\lambda_{j}|\leq\sum_{\ell=0}^{N-1}G_{\ell}$ for all $j$. Since

[TABLE]

we have that

[TABLE]

The strict inequality is due to $-1\leq erf(x)\leq 1$ for all $x$ and $N$ being finite. It follows that the eigenvalues of $A^{T}A$ are $|\lambda_{j}|^{2}<1$. Therefore, the spectrum of $A^{T}A-I$ is {\color[rgb]{0,0,0}s\normalcolor}(A^{T}A-I)<0. ∎

Therefore, $E_{k}^{m+1}-E_{k}^{m}\leq 0$. Consequently,

[TABLE]

which proves that the numerical method is stable.

6 Numerical results

6.1 The 1-dimensional model

For our 1-dimensional simulations, we perform 100,000 realizations of the ABM with agent $\mathfrak{u}$ initialized at $x_{0}=0.5$. The agent moves with spatial step size of $\triangle x=0.01$ and time step $\triangle t=\triangle x^{2}/2$. For the corresponding PDE model, we use the stable numerical algorithm detailed in Section 5 with a point source at $x_{0}=0.5$. We choose $N$ so that $G_{0},G_{N-1}<\varepsilon_{mach}$ and we assign the numerical step sizes as $dx=\triangle x$, $dt=\triangle t$, and $d\xi=\xi_{c}/2000$. In both the ABM and PDE model, we define the agent absorption function as $\beta(x)=\alpha\int_{B(x,\triangle x/2)}C(x)\,dx$. The $\alpha$ parameter defines the permeability of the agent’s membrane and for the following examples, we let $\alpha=0.1$.

6.1.1 Example 1: Max concentration at starting location

For this example, the chemical concentration $C(x)=\frac{1}{1+10(x-0.5)^{2}}$ is symmetric and concave down around $x=0.5$. A comparison of the ABM and our continuum PDE model is shown in Fig. 2, for a critical or tolerance threshold of $\xi_{c}=10\triangle x\triangle t$. The distribution of cells or agents in the initial live state are shown in color with time on the vertical axis and spatial location on the horizontal axis. The values on Fig. 2b at location $(x_{i},t_{m})$ are the numerical solutions $p_{i}^{m}$ from (30), which are interpreted as the probability a cell is alive and located within region $B(x_{i},\triangle x/2)$ at time $t$. Since the agents are all initialized at $x_{o}=0.5$, we observe a high density of cells close to this point for small time intervals. We note that Fig. 2b is smoother than 2a since it is a continuous approximation whereas the ABM has agents moving discretely either to the left or right at each time step.

Additionally, since $C(x)$ has a max at $x=0.5$, this causes the probability distribution $p(x,t)$ to become bimodal at approximately $t=0.0055$. Those cells that have remained close to the initial starting location have absorbed more particles than those that have moved left or right. Hence, cells close to $x=0.5$ are moving out of the initial live cell state when they reach their absorption capacitance $\xi_{c}$.

The probability an agent is alive at a given time $t$ is the survival probability $P(t)$, calculated as

[TABLE]

In Fig. 3a, we observe that $P(t)$ for the ABM simulation and PDE approximations match; there is a sharp decrease in survival probability after $t=0.005$ and the majority of the cells have died at $t=0.007$.

The mean location of the live agents is calculated as $\mu(t)=\int_{\mathbb{R}}x\hat{p}(x,t)\,dx$, where $\hat{p}(x,t)=p(x,t)/P(t)$ is the normalized value of $p(x,t)$ at each time $t$. The numerical PDE solution solves for the average value in the interval centered at $x_{i}$ with radius $\triangle x/2$, $B(x_{i},\triangle x/2)$. This allows the calculation of $\mu(t)$, the mean at time $t=m\tau$, as

[TABLE]

the exact integral of the approximate piece-wise constant solution. Just as we did when calculating the convolution, we can take $p_{i}^{m}$ out of the integral since it is piece-wise constant. In a similar way, we can calculate $\sigma^{2}(t)$, the variance at time $t=m\tau$, as

[TABLE]

The mean location of the ABM simulation and PDE approximation is shown in Fig. 3b. The chemical concentration $C(x)$ is symmetric around $x=0.5$, the location where the agents are initialized, and there is no bias in movement ($\ell(x)=r(x)=0.5$). Hence, we would expect the mean location of agents in the initial state to be centered at $x=0.5$. We see that until approximately $t=0.006$, the PDE mean and the ABM mean are close to $x=0.5$. For times $t>0.006$, the number of agents in the ABM simulation is relatively small, as shown in Fig. 3a. This accounts for the increasing stochastic noise in the mean, as well as the standard deviation, which is shown in Fig. 3c.

At each iteration of the ABM simulation, the agent can move either left or right. We see that the agents that remain in the initial state are those that are furthest from $x=0.5$, where $C(x)$ is smaller than at $x=0.5$. As a result, the standard deviation is a monotonically increasing function, as seen in Fig. 3c. At approximately $t=0.005$, many cells towards the center of the simulation change state, which causes the “corner” in the standard deviation graph.

6.1.2 Example 2: Decreasing concentration

The chemical concentration is $C(x)=\exp\left(-x^{2}\right)$, which is monotonically decreasing in the interval $[0,1]$ and all agents or cells are initialized at $x_{o}=0.5$. We expect that the agents which tend to move to the right within this interval have a higher probability of remaining in the initial state. As shown in Fig. 4, the cells that remain in the initial state tend to be further to the right and again, we have excellent qualitative agreement between the ABM and the new PDE continuum model. In Fig. 4a we observe a striped pattern, which is a result of the ABM agents moving only left or right at any given iteration. At a critical threshold of $\xi_{c}=10\triangle x\triangle t$, cells are able to achieve a cumulative chemical absorption $\xi>\xi_{c}$, causing the cell to transition states or die. The survival probability shows this trend in Fig. 5a, where there is a sharp decrease in survival probability at $t=0.055$.

To further characterize the agreement between the ABM simulation and our PDE approximation, we again look at the mean and standard deviation of the location of live cells (with cumulative absorption $\xi<\xi_{c}$). In Fig. 5b, we observe that the mean location (calculated using (35)) does move to the right of the initial location $x_{o}=0.5$ due to the decreased concentration $C(x)$ to the right of $x=0.5$ (allowing cells to live in this region for a longer period of time). Again, we see that there is noise in the ABM mean for times $t>0.008$, when there are relatively few agents in the initial state.

As shown in Fig. 5c, the standard deviation of the agents locations is increasing for $0\leq t\leq 0.005$, which corresponds to the time interval where most cells are alive (see survival probability in Fig. 5a). At $t=0.005$, agents with a cumulative absorption reaching $\xi_{c}$ begin to change state. Cells to the right of $x_{o}=0.5$ tend to remain in the initialized state, which moves the mean to the right and reduces the variance. A majority of the cells have changed state by $t=0.008$, where the cells that remain are those that continued to move right. Thus, the standard deviation approaches zero. Similar to Example 1, we see that as the number of agents in the ABM simulation approaches zero, the stochastic noise influences the variance (Fig. 5c).

6.1.3 Example 3: Biased random walk

Suppose the random walk has a constant bias, where $\ell$ and $r$ denote the probabilities of moving left or right, respectively. Our absorption model is the following PDE

[TABLE]

where $a=\triangle x(\ell-r)/\triangle t$ and $D=\triangle x^{2}(\ell+r)/(2\triangle t)$. Note that the existence proof also holds if the agent moves with a constant bias. After splitting the linear operator, the only difference between (37) and our initial cumulative absorption model equation (5) is the form of the Green’s function

[TABLE]

Replacing the diffusion Green’s function with this advection-diffusion Green’s function does not affect the existence proof in Section 4.3. Further, by integration by parts and using our far-field boundary condition, we can show that $\int_{\Omega}aUU_{x}\,d\xi\,dx=0$. Therefore, the biased model (37) satisfies Theorem 1, so we can prove that the PDE (37) is well-posed.

We set the chemical concentration as $C(x)=1/(1+10(x-0.5)^{2})$ and the absorption capacitance as $\xi_{c}=10\triangle x\triangle t$, the same as in Example 1. However, in contrast to Example 1, we set the probability an agent moves left as $\ell=0.4$ and the probability an agent moves right as $r=0.6$. We note that a larger density of agents tend to move to the right. In fact, the graph in Fig. 7b initially moves to the right in a straight line at a rate of $0.2\triangle x$. This is due to the fact that only the biased movement determines the agent locations. The graph of the mean location makes a sudden change around $t=0.005$, which is when some agents absorb above the threshold $\xi_{c}$ and begin leaving the live state. At that time, the chemical profile begins to influence the mean location of agents, which also explains the shapes of the distributions in Fig. 6.

The survival probability for the ABM and continuum PDE model is shown in Fig. 7a, and again there is good agreement between the ABM and PDE solution. In comparison to the unbiased movement case in Fig. 3a, Fig. 7a with biased movement begins decreasing at an earlier time, but then decreases at a slower rate.

6.2 The 2-dimensional model

We can readily extend the analysis and numerical methods in Sections 4-5 to the 2-dimensional case. To account for the increased stochasticity of adding an additional dimension, we initialize 10 million agents. The agents in the ABM move with spatial step size of $\triangle x=\triangle y=0.01$ and time step $\triangle t=\triangle x^{2}/2$. Similarly, the PDE model utilizes a spatial step size of {\color[rgb]{0,0,0}dx=dy=\triangle x\normalcolor} and a time step of $dt=\triangle t$, and cumulative absorption of $d\xi=\xi_{c}/1000$. For both the ABM and PDE model, we set $\beta(\mathbf{x})=\alpha\int_{B(\mathbf{x},\triangle x/2)}C(\mathbf{x})\,d\mathbf{x}$ where the chemical concentration is $C(x,y)=0.5(\sin(4\pi x)\sin(4\pi y)+1)$ and the chemical absorption threshold is $\xi_{c}=2\triangle x\triangle y\triangle t$.

The surface plot of the concentration local to the initialized agents in $[0,1]\times[0,1]$ is shown in the dashed line contour plots in Fig. 8 and Fig. 9, where lighter colored lines denote a value closer to 0 and darker colored lines denote a value closer to 1. The concentration is symmetric along the lines $y=x$ and $y=1-x$. Near the initial location at $(0.5,0.5)$, there are local concentration minimums along the line $y=-x$. Thus, it makes sense that the probabilities for agents in the initial live state tend to be higher close to these chemical sinks, as shown in Figs. 8 and 9. In fact, Figs. 8(b)-(c) and 9(b)-(c) show the probability density function mode bifurcation. That is, the chemical distribution causes $p_{i}^{m}$ to evolve into a bi-modal distribution, with each peak located on the line $y=1-x$ and equidistant to the line $y=x$. Again, when comparing the survival probability as a function of time, we observe excellent agreement between the ABM and continuum PDE (Fig. 10a).

Fig. 10b demonstrates that the mean location of the PDE approximation remains constant at $(0.5,0.5)$. The ABM mean is not constant. However, since the ABM mean is contained within the region $B((0.5,0.5),\triangle x/2)$, and travels away from the PDE mean for times $t>0.003$, we can assume that this is due to the greater influence of stochastic noise as the number of agents in the initial state becomes relatively small. Since there are sufficiently many agents towards the end of the simulation and the mean during this simulation is within the control region $B((0.5,0.5),\triangle x/2)$, we see in Fig. 10c the standard deviation of the ABM data is not unduly influenced by the stochastic noise. Hence, the ABM and PDE standard deviation curves match reasonably well throughout the simulations.

We can see the difference in how the model develops if we decrease the absorption proportion parameter to $\alpha=0.01$ in Fig. 11 and Fig. 12. The agents diffuse for a longer time period before absorbing sufficient chemical to split. With $\alpha=0.01$, the distribution forms two peaks around $t=0.02$ (as shown in Fig. 11(d), 12(c)), whereas with $\alpha=0.10$, the distribution forms two peaks around $t=0.0024$ (as shown in Fig. 8(c), 9(c)). The pattern is different from $\alpha=0.10$ as $t$ increases further in that since $p_{i}^{m}$ initially diffuses farther before changing states, the distribution will settle along additional chemical sinks (as shown in Fig. 11(d), 12(d)).

If we instead choose an $\alpha\ll 1$, the diffusion time of the agents is much faster than the time for the agents to absorb chemical to capacitance. In that case, we can simplify and re-frame the model as a diffusion-dominant absorption model. The derivation and examples of such a model are further developed in Section 7.

7 Relative scaling approximations

The above examples assumed a particular scaling of parameters in (5) and (6). However, other scaling relations may be possible and we can further simplify the equations or analysis. First, let us non-dimensionalize the absorption variable, $\xi$, by using the scaling factor $\xi_{c}$ to obtain $\bar{\xi}=\xi/\xi_{c}$. We can rewrite (5) as

[TABLE]

with the agent changing state when $\xi\geq 1$. There are three different regimes based on the term $\beta(x)/(D\xi_{c})$.

First, we will investigate when $\beta(x)/(D\xi_{c})\gg 1$. This occurs when the agent absorbs chemical above its capacitance before diffusion moves the agent, as can be observed in Fig 13(a). Assuming $\beta(x)>0$ is continuous, we can asymptotically simplify (38) to

[TABLE]

with a solution $U(x,t,\bar{\xi})=\phi(x,\bar{\xi}-\beta(x)t/\xi_{c})$ for all $x\in\mathbb{R}$. An example solution of the total density in the live state in an absorption-dominant parameter regime is shown in Fig 13(b). Note that as $\langle\beta(x)\rangle/(D\xi_{c})$ approaches $\infty$, the ABM and PDE densities converge.

Second, we will investigate when $\beta(x)/(D\xi_{c})\ll 1$. This occurs when the agent diffuses much faster than the agent absorbs chemical, as is shown in Fig 14(a). We cannot make an asymptotic argument, ignoring the absorption term, since that is our primary interest in this model. Thus, we make the assumption that the density of the agents in the live state quickly becomes uniform over the spatial coordinate. This assumption allows us to collapse the spatial coordinate and have a solution solely in absorption and time dimensions. We can define the PDE initial condition as $\bar{\phi}(\bar{\xi}):=\int_{\mathbb{R}}\phi(x,\bar{\xi})\,dx$ and we can redefine the PDE absorption term as the average value of $\beta$ in the spatial domain, $\langle\beta(x)\rangle$, to obtain the following PDE,

[TABLE]

The solution is $U(\bar{\xi},t)=\bar{\phi}\left(\bar{\xi}-\langle\beta(x)\rangle t/\xi_{c}\right)$. An example solution of the total density in the live state in a diffusion-dominant parameter regime is shown in Fig 14(b). Note that as $\langle\beta(x)\rangle/(D\xi_{c})$ approaches 0, the ABM and PDE densities converge.

8 Discussion and Conclusions

In this work, we have developed a continuum PDE approximation to a stochastic agent-based model that has a cumulative coupling to the environment. We have shown through simulations that the ABM agrees qualitatively with the governing PDE. We have analyzed the newly developed PDE, showing that we have developed a stable, well-posed equation.

The modeling framework developed assumes that the cells have cumulative exposure but the chemical or substance is not being removed from the environment. This setup can be used to model morphogens (signaling molecules) that often act directly on a cell by binding to a receptor. It is well established that cells exposed to high levels of transforming growth factor (TGF)-$\beta$ will lead to cell death [10]. In this scenario, morphogens bind and initiate a secondary process in the cell that accumulates and leads to cell death. After a small period of time (smaller than the scale of movement or time to cell death), the morphogen is released from the receptor and thus the relative morphogen concentration can be assumed constant in time.

Now that we have developed the initial framework, additional model extensions will allow us to look at real world examples for cells coupled to the environment in a cumulative way. There are many examples of cells that have directed motion due to a chemoattractant [28]. This would involve extending the model beyond the general framework that studied a constant or fixed bias without affecting the chemical concentration (Section 6.1.3). Cell movement will be biased due to the chemical gradient and, in turn, the chemical gradient will change due to localized cellular absorption. We can derive the system of PDEs in a manner similar to that of Sections 3.1 and 3.2. The method of deriving the Heaviside term for $V$ in Section 3.2 provides us the capability of handling different states as well as the capability of handling cells that die and no longer absorb chemicals. However, further analytical work will be required since the well-posedness proof of the PDE and numerical solution will need to be modified to account for a new space and time dependent advection term as well as the fact that the absorption term, $\beta$, is coupled with the space and time dependent chemical profile.

Further, there are other interesting biological examples of toxins that might not only kill the cell, but will change their motility. In wound healing, the migration of epithelial cells is important. However, exposure to chemical substances (e.g. chemical warfare agents) has been shown to also decrease cell migration speeds [27]. Another example is in cancer therapies where immunotoxins such as MOC31PE reduce cell migration and induce gene expression, leading to cell death in ovarian cancer cells [32]. Additionally, human neural crest cells migrate to particular locations in the embryo and differentiate into different cell types. It has been found that interferon-beta (IFN$\beta$), a potential toxin [21], impairs neural crest cell migration. To model these types of applications, one would potentially have to re-derive the equations for a variable step size based on the cumulative exposure. In terms of the computational method, we believe this case could be simulated using the same semi-discrete numerical formulation since the $\xi$-dependent diffusion parameter will be treated as a constant in the Green’s function. However, the $\xi$-dependent diffusion parameter will cause subtle complications when manipulating the integrals within the well-posedness analysis.

Although we have focused on the example of a cell or agent absorbing a fraction of the particles of a chemical in the surrounding environment, these equations are generally applicable to any scenario where a cell is accumulating any quantity that is a function of space and/or time. In the current contexts, we have focused on the case of cells that do not interact with each other. The focus of future work will investigate interactions of different cells or agents as they interact in a cumulative way with their environment.

Acknowledgments

This work was supported, in part, by National Science Foundation Division of Mathematical Sciences grant 1455270. S.D. Olson was also supported, in part, by a Fulbright research scholar grant. M.A. Yereniuk has been supported, in part, by the SMART scholarship program. All simulations were run on a cluster acquired through National Science Foundation grant 1337943.