Reactions, Diffusion and Volume Exclusion in a Heterogeneous System of Interacting Particles

TL;DR

This paper introduces a new reaction-diffusion modeling framework that incorporates volume exclusion effects for heterogeneous particles, bridging microscopic interactions and macroscopic descriptions, with applications to chemotaxis.

Contribution

The authors develop a novel macroscopic PDE system derived from microscopic models that explicitly accounts for volume exclusion and reactions among heterogeneous particles.

Findings

Reaction terms are of lower order than nonlinear diffusion in the asymptotic expansion.

Including the next reaction term improves agreement with microscopic simulations.

The macroscopic model facilitates direct parameterization from experimental data.

Abstract

Complex biological and physical transport processes are often described through systems of interacting particles. Excluded-volume effects on these transport processes are well studied, however the interplay between volume exclusion and reactions between heterogenous particles is less well known. In this paper we develop a novel framework for modeling reaction-diffusion processes which directly incorporates volume exclusion. From an off-lattice microscopic individual based model we use the Fokker--Planck equation and the method of matched asymptotic expansions to derive a low-dimensional macroscopic system of nonlinear partial differential equations describing the evolution of the particles. A biologically motivated, hybrid model of chemotaxis with volume exclusion is explored, where reactions occur at rates dependent upon the chemotactic environment. Further, we show that for reactions due to contact interactions the appropriate reaction term in the macroscopic model is of lower order in the asymptotic expansion than the nonlinear diffusion term. However, we find that the next reaction term in the expansion is needed to ensure good agreement with simulations of the microscopic model. Our macroscopic model allows for more direct parameterization to experimental data than the models available to date.

Equations

$$\mathrm{d}{\mathbf{X}}_i(t)$$

$${\partial }_t\left(\begin{array}{c}b\\ r\end{array}\right)$$

$$\mathcal{D}=\text{diag}(D_{\mathsf{b}},D_{\mathsf{r}})+{ϵ}^d\tilde{\mathcal{D}}(b,r),$$

$$\tilde{\mathcal{D}}=\left(\begin{array}{c}{D}_{\mathsf{b}}(N_{\mathsf{b}}-1)\omega b-D_{\mathsf{b}}{N}_{\mathsf{r}}{\eta }_{\mathsf{b}}r{D}_{\mathsf{b}}{N}_{\mathsf{r}}{\xi }_{\mathsf{b}}b\\ \\ D_{\mathsf{r}}{N}_{\mathsf{b}}{\xi }_{\mathsf{r}}r{D}_{\mathsf{r}}(N_{\mathsf{r}}-1)\omega r-D_{\mathsf{r}}{N}_{\mathsf{b}}{\eta }_{\mathsf{r}}b\end{array}\right),$$

$$\mathcal{F}=\left(\begin{array}{cc}{\mathbf{f}}_{\mathsf{b}}& N_{\mathsf{r}}{ϵ}^d{\eta }_{\mathsf{b}}({\mathbf{f}}_{\mathsf{r}}-{\mathbf{f}}_{\mathsf{b}})b\\ N_{\mathsf{b}}{ϵ}^d{\eta }_{\mathsf{r}}({\mathbf{f}}_{\mathsf{b}}-{\mathbf{f}}_{\mathsf{r}})r& {\mathbf{f}}_{\mathsf{r}}\end{array}\right),$$

$$\begin{array}{rl}\omega & =\frac{2\pi (d-1)}{d},{\xi }_k=\frac{2\pi }{d}\frac{[(d-1)D_k+d{D}_l]}{D_k+D_l},\\ {\eta }_k& =\frac{2\pi }{d}\frac{D_k}{D_k+D_l}.\end{array}$$

$$\mathrm{d}{\mathbf{X}}_i(t)$$

$$B\stackrel{k_{\mathsf{b}}}{\to }R,R\stackrel{k_{\mathsf{r}}}{\to }B,$$

$${\mathcal{B}}_{ϵ}=\left\{\vec{x}\in {\Omega }^N:\exists i\ne j\text{ such that }\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert \le ϵ\right\}.$$

$$\begin{array}{rl}{\partial }_{t}P=& {\nabla }_{\vec{x}}\cdot \left[{\mathsf{D}}_{\vec{s}}{\nabla }_{\vec{x}}P(\vec{x},\vec{s},t)-\vec{F}(\vec{x},\vec{s})P(\vec{x},\vec{s},t)\right]\\ & +k_{\mathsf{b}}\sum _{\vec{c}\in \mathcal{B}(\vec{s})}P(\vec{x},\vec{c},t)+k_{\mathsf{r}}\sum _{\vec{c}\in \mathcal{R}(\vec{s})}P(\vec{x},\vec{c},t)\\ & -\left(k_{\mathsf{b}}|\mathcal{R}(\vec{s})|+k_{\mathsf{r}}|\mathcal{B}(\vec{s})|\right)P(\vec{x},\vec{s},t),\end{array}$$

$$\left[{\mathsf{D}}_{\vec{s}}{\nabla }_{\vec{x}}P(\vec{x},\vec{s},t)-\vec{F}(\vec{x},\vec{s})P(\vec{x},\vec{s},t)\right]\cdot \vec{n}=0,$$

$$b(\mathbf{x},t)={\int }_{{\Omega }_{ϵ}^N}\sum _{\vec{s}\in {\mathcal{C}}^N}P(\vec{x},\vec{s},t)\delta (s_1-\mathsf{b})\delta (\mathbf{x}-{\mathbf{x}}_1)\mathrm{d}\vec{x},$$

$$r(\mathbf{x},t)={\int }_{{\Omega }_{ϵ}^N}\sum _{\vec{s}\in {\mathcal{C}}^N}P(\vec{x},\vec{s},t)\delta (s_1-\mathsf{r})\delta (\mathbf{x}-{\mathbf{x}}_1)\mathrm{d}\vec{x}.$$

$${\partial }_t{P}_{\mathsf{bb}}$$

$${\partial }_t{P}_{\mathsf{br}}$$

$${\partial }_t{P}_{\mathsf{rb}}$$

$${\partial }_t{P}_{\mathsf{rr}}$$

$$b({\mathbf{x}}_1,t)={\int }_{{\Omega }_{ϵ}({\mathbf{x}}_1)}(P_{\mathsf{bb}}+P_{\mathsf{br}})\mathrm{d}{\mathbf{x}}_2,$$

$$r({\mathbf{x}}_1,t)={\int }_{{\Omega }_{ϵ}({\mathbf{x}}_1)}(P_{\mathsf{rr}}+P_{\mathsf{rb}})\mathrm{d}{\mathbf{x}}_2,$$

$$\begin{array}{rl}{\partial }_{t}b=& {\int }_{{\Omega }_{ϵ}({\mathbf{x}}_1)}{\nabla }_{{\mathbf{x}}_1}\cdot \left(D_{\mathsf{b}}{\nabla }_{{\mathbf{x}}_1}{P}_{\mathsf{bb}}+D_{\mathsf{b}}{\nabla }_{{\mathbf{x}}_1}{P}_{\mathsf{br}}\right)\mathrm{d}{\mathbf{x}}_2\\ & +{\int }_{{\Omega }_{ϵ}({\mathbf{x}}_1)}{\nabla }_{{\mathbf{x}}_2}\cdot \left(D_{\mathsf{b}}{\nabla }_{{\mathbf{x}}_2}{P}_{\mathsf{bb}}+D_{\mathsf{r}}{\nabla }_{{\mathbf{x}}_2}{P}_{\mathsf{br}}\right)\mathrm{d}{\mathbf{x}}_2\\ & +k_{\mathsf{r}}r-k_{\mathsf{b}}b.\end{array}$$

$$\displaystyle\begin{aligned} \partial_{t}b&=D_{\mathsf{b}}\nabla_{{\bf x}_{1}}^{2}b+k_{\mathsf{r}}r-k_{\mathsf{b}}b\\ &\ +\!\int_{\partial B_{\epsilon}({\bf x}_{1})}\!-\big{[}2D_{\mathsf{b}}\nabla_{{\bf x}_{1}}(P_{\mathsf{b}\mathsf{b}}+P_{\mathsf{b}\mathsf{r}})\\ &\qquad\qquad\qquad+(D_{\mathsf{b}}-D_{\mathsf{r}})\nabla_{{\bf x}_{2}}P_{\mathsf{b}\mathsf{r}}\big{]}\cdot\mathbf{n_{2}}\mathrm{d}S_{{\bf x}_{2}},\end{aligned}$$

$$\displaystyle\begin{aligned} \partial_{t}r&=D_{\mathsf{r}}\nabla_{{\bf x}_{1}}^{2}r-k_{\mathsf{r}}r+k_{\mathsf{b}}b\\ &\ +\!\int_{\partial B_{\epsilon}({\bf x}_{1})}\!-\big{[}2D_{\mathsf{r}}\nabla_{{\bf x}_{1}}(P_{\mathsf{r}\mathsf{b}}+P_{\mathsf{r}\mathsf{r}})\\ &\qquad\qquad\qquad+(D_{\mathsf{r}}-D_{\mathsf{b}})\nabla_{{\bf x}_{2}}P_{\mathsf{r}\mathsf{b}}\big{]}\cdot\mathbf{n_{2}}\mathrm{d}S_{{\bf x}_{2}}.\end{aligned}$$

$$P_{\mathsf{bb}}^{\text{out}}=q_{\mathsf{b}}({\mathbf{x}}_1,t)q_{\mathsf{b}}({\mathbf{x}}_2,t)+ϵP_{\mathsf{bb}}^{\text{out},(1)}({\mathbf{x}}_1,{\mathbf{x}}_2,t)+\cdots ,$$

$$\displaystyle\epsilon^{2}\partial_{t}\hat{P}_{\mathsf{b}\mathsf{b}}=\$$

$$\displaystyle+\epsilon^{2}\big{[}D_{\mathsf{b}}\nabla_{\hat{{\bf x}}_{1}}^{2}\hat{P}_{\mathsf{b}\mathsf{b}}+k_{\mathsf{r}}(\hat{P}_{\mathsf{b}\mathsf{r}}+\hat{P}_{\mathsf{r}\mathsf{b}})-2k_{\mathsf{b}}\hat{P}_{\mathsf{b}\mathsf{b}}\big{]}.$$

$$2\hat{\mathbf{x}}\cdot {\nabla }_{\hat{\mathbf{x}}}{\hat{P}}_{\mathsf{bb}}=ϵ\hat{\mathbf{x}}\cdot {\nabla }_{{\hat{\mathbf{x}}}_1}{\hat{P}}_{\mathsf{bb}},\text{on}\Vert \hat{\mathbf{x}}\Vert =1.$$

$$\displaystyle\begin{aligned} \hat{P}_{\mathsf{b}\mathsf{b}}\sim\ &q_{\mathsf{b}}^{2}(\hat{{\bf x}}_{1},t)+\epsilon\big{[}q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)\hat{{\bf x}}\cdot\nabla_{\hat{{\bf x}}_{1}}q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)\\ &+P_{\mathsf{b}\mathsf{b}}^{\text{out},(1)}(\hat{{\bf x}}_{1},\hat{{\bf x}}_{1},t)\big{]},\end{aligned}$$

$${\nabla }_{\hat{\mathbf{x}}}^2{\hat{P}}_{\mathsf{bb}}^{(0)}$$

$$\hat{\mathbf{x}}\cdot {\nabla }_{\hat{\mathbf{x}}}{\hat{P}}_{\mathsf{bb}}^{(0)}$$

$${\hat{P}}_{\mathsf{bb}}^{(0)}$$

$${\nabla }_{\hat{\mathbf{x}}}^2{\hat{P}}_{\mathsf{bb}}^{(1)}$$

$$\hat{\mathbf{x}}\cdot {\nabla }_{\hat{\mathbf{x}}}{\hat{P}}_{\mathsf{bb}}^{(1)}$$

$${\hat{P}}_{\mathsf{bb}}^{(1)}$$

$$+P_{\mathsf{bb}}^{\text{out},(1)}({\hat{\mathbf{x}}}_1,{\hat{\mathbf{x}}}_1,t),\text{as}\Vert \hat{\mathbf{x}}\Vert \to \infty ,$$

$$\displaystyle\begin{aligned} \hat{P}_{\mathsf{b}\mathsf{b}}\sim\ &q_{\mathsf{b}}^{2}(\hat{{\bf x}}_{1},t)+\epsilon\big{[}q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)\hat{{\bf x}}\cdot\nabla_{\hat{{\bf x}}_{1}}q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)\\ &+P_{\mathsf{b}\mathsf{b}}^{\text{out},(1)}(\hat{{\bf x}}_{1},\hat{{\bf x}}_{1},t)\big{]}.\end{aligned}$$

$$\displaystyle\begin{aligned} \hat{P}_{\mathsf{r}\mathsf{r}}\sim q_{\mathsf{r}}^{2}(\hat{{\bf x}}_{1},t)+\epsilon\big{[}&q_{\mathsf{r}}(\hat{{\bf x}}_{1},t)\hat{{\bf x}}\cdot\nabla_{\hat{{\bf x}}_{1}}q_{\mathsf{r}}(\hat{{\bf x}}_{1},t)\\ &+P_{\mathsf{r}\mathsf{r}}^{\text{out},(1)}(\hat{{\bf x}}_{1},\hat{{\bf x}}_{1},t)\big{]}.\end{aligned}$$

$$P_{\mathsf{br}}^{\text{out}}=q_{\mathsf{b}}({\mathbf{x}}_1,t)q_{\mathsf{r}}({\mathbf{x}}_2,t)+ϵP_{\mathsf{br}}^{\text{out},(1)}({\mathbf{x}}_1,{\mathbf{x}}_2,t)+\cdots ,$$

$${\hat{P}}_{\mathsf{br}}({\hat{\mathbf{x}}}_1,\hat{\mathbf{x}},t)=P_{\mathsf{br}}({\mathbf{x}}_1,{\mathbf{x}}_2,t)$$

$${ϵ}^2{\partial }_t{\hat{P}}_{\mathsf{br}}$$

$$\displaystyle+\epsilon^{2}\big{[}D_{\mathsf{b}}\nabla_{\hat{{\bf x}}_{1}}^{2}\hat{P}_{\mathsf{b}\mathsf{r}}+k_{\mathsf{b}}\hat{P}_{\mathsf{b}\mathsf{b}}+k_{\mathsf{r}}\hat{P}_{\mathsf{r}\mathsf{r}}-(k_{\mathsf{b}}\!+\!k_{\mathsf{r}})\hat{P}_{\mathsf{b}\mathsf{r}}\big{]}.$$

Full text

I Introduction

Cellular migration plays an important role in many biological processes, including tumor growth and invasion through an extracellular matrix Chauviere et al. (2010); Corcoran and Maestro (2003), the formation of blood vessels via the movement of endothelial cells during embryogenesis Drake (2003); Risau (1997) and the directed motion of immune cells to infected sites Robertson et al. (1993); Wyckoff et al. (2004). Classical continuum models assume that cells move down spatial gradients in their density (diffusion) and/or up spatial gradients in the concentration of a diffusible chemical (chemotaxis). However, these models treat the cells as point particles. In recent years experimental evidence has revealed the importance of excluded-volume effects in a variety of biological applications such as cellular migration Mort et al. (2016); Kim and Yethiraj (2010); Gejji et al. (2012); Plank and Simpson (2012) and molecular traffic Neri et al. (2011); Graf and Frey (2017); Chowdhury et al. (2005), as well as physical processes including vehicular traffic Chowdhury et al. (2000) and queueing Kipnis (1986). These findings have motivated a significant amount of theoretical research focussed on incorporating excluded-volume effects into mathematical models of transport Sun and Weinstein (2007); Lushnikov et al. (2008); Pillay et al. (2018).

One approach is to discretize space into a regular lattice, and allow at most one particle to occupy each lattice site Simpson et al. (2009); Ascolani et al. (2013). Diffusion is represented by jumps between neighboring lattices and these are aborted if the neighbouring lattice site is already occupied. An alternative lattice-based approach which incorporates finite-size effects is the cellular Potts model Lushnikov et al. (2008). An advantage of these lattice-based approaches is that it is relatively straightforward to include complex individual-based mechanisms, and to obtain a macroscopic model based on partial differential equations (PDEs). However, the form of the lattice can introduce artefacts into these models, such as a bias towards homogeneity Flache and Hegselmann (2001); Penington et al. (2011).

An alternative approach is to derive continuum equations from an off-lattice individual-based model, starting from Brownian particles with excluded-volume interactions. For a single species of particles one obtains a nonlinear diffusion equation Bruna and Chapman (2012a), whereas on-lattice models yield a linear diffusion equation in the continuum limit Simpson et al. (2009). An off-lattice model with volume exclusion, but with motility rules inherited from the standard lattice-based approach is considered in Dyson et al. (2012). In this case, the resulting PDE model is also a nonlinear diffusion equation similar to Bruna and Chapman (2012a).

As in the case of diffusion only, stochastic reaction-diffusion processes have been modeled using on- and off-lattice approaches Meinecke and Eriksson (2016); Seki et al. (2012). A challenge with lattice-based approaches is that it is difficult to combine reactions and finite-size interactions. As mentioned above, in the latter case the lattice size is chosen so that at most one cell occupies a single lattice site. By contrast, the usual lattice-based approach for bimolecular reactions requires two cells to be in the same compartment in order to react Erban and Chapman (2009). The standard off-lattice approach to model reaction-diffusion processes is based on the Smoluchowski model, which represents cells as point particles undergoing Brownian motion and reacting whenever they become closer than a given reaction radius Andrews and Bray (2004). The reaction radius is chosen as a function of the reaction rate and diffusion constants of the reactants, and is typically independent of their physical radii. In fact, it has been shown that the reaction radius may be unrealistically smaller than the radii of the reactants Erban and Chapman (2009). Several methods have been proposed to overcome the shortcomings of these models Andrews and Bray (2004); Isaacson (2013), including one where the reaction radius can be fixed to represent the reactant radius Erban and Chapman (2009), but to our knowledge none includes the excluded-volume effects which have been shown to significantly impact diffusion-limited reactions Park et al. (2003). To this end, in this paper we extend the off-lattice model in Bruna and Chapman (2012a) to account for reactions between two subpopulations of hard-sphere Brownian particles. The advantage of our approach to modelling reactions is that the sizes of the particles (and therefore excluded-volume interactions) are included explicitly in our individual-based model, thereby circumventing the problem of unrealistically small reaction radii and providing a more natural way to describe diffusion-limited reactions Erban and Chapman (2009).

The remainder of the paper is organised as follows. In Section II we introduce unimolecular reactions to represent spontaneous switching between subpopulations. These simple reactions allow us to develop our framework, but can also be motivated by the phenotypic switching of tumor cells within a glioblastoma Corcoran and Maestro (2003) as well as other canonical studies in mathematical biology such as infectious diseases and social interaction networks Murray (2002). We first formulate the discrete model for the interacting particles and then use the Fokker–Planck equation in combination with the method of matched asymptotic expansions to derive a low-dimensional system of nonlinear reaction-diffusion PDEs. A numerical example comparing the microscopic and macroscopic models is presented in Subsection II.6. In Subsection II.7 we apply our framework to cellular chemotaxis, by making the rate at which cells switch between phenotypes depend on the concentration of a third species which also acts as a chemotattractant for one cell population. In Section III we study the case when reactions between subpopulations (bimolecular reactions) occur in response to collisions. A biological motivation for switching due to local interactions can be found in juxtacrine signaling as illustrated by the delta-notch signaling pathway Androutsellis-Theotokis et al. (2006). Finally in Section IV we summarize our results and suggest possible directions for future investigations.

II Spontaneous Switching: Unimolecular Reactions

We place our work in context by briefly summarizing the results of Bruna and Chapman (2012b), which considers two subpopulations of inert Brownian particles with hard-core excluded-volume interactions. Consider a population of $N$ hard-spheres of diameter $\epsilon\ll 1$ within a bounded domain $\Omega\subset\mathbb{R}^{d}$ of typical dimensionless volume of order one, where $d=2,3$. There are two subpopulations of spheres, which we term “red” and “blue”. Let there be $N_{\mathsf{b}}$ blue spheres and $N_{\mathsf{r}}$ red spheres, where $N=N_{\mathsf{b}}+N_{\mathsf{r}}$. Assume that the particles occupy a small volume fraction, so that $(N_{\mathsf{b}}+N_{\mathsf{r}})\epsilon^{d}\ll 1$. Let $\mathcal{I}_{\mathsf{b}}$ denote the set of indices corresponding to blue particles, and $\mathcal{I}_{\mathsf{r}}$ the indices corresponding to red particles. The blue and red spheres have diffusion coefficients $D_{\mathsf{b}}$ and $D_{\mathsf{r}}$, and drift vectors ${\mathbf{f}}_{\mathsf{b}}$ and ${\mathbf{f}}_{\mathsf{r}}$, respectively. The center of the $i$th particle is denoted by ${\mathbf{X}}_{i}(t)$ and evolves according to the overdamped Langevin SDE:

[TABLE]

for $1\leq i\leq N$, where the ${\mathbf{W}}_{i}$ are independent, $d$-dimensional standard Brownian motions. Reflective boundary conditions are imposed whenever two spheres are in contact ($\|{\bf X}_{i}-{\bf X}_{j}\|=\epsilon$ for $i\neq j$) and also on the boundary of the domain $\partial\Omega$. From a particle-level description, the method of matched asymptotic expansions for a small but finite volume fraction is used to obtain a continuum model for the marginal probability density functions, $b({\mathbf{x}},t)$ and $r({\mathbf{x}},t)$ for blue and red spheres, respectively. The continuum model comprises the following system of nonlinear cross-diffusion equations

[TABLE]

where the constants $\omega,\xi_{k},\eta_{k}$ (for $k=\mathsf{b},l=\mathsf{r}$ and vice versa) are defined by

[TABLE]

The density-dependent terms in (2b), contained in $\widetilde{\mathcal{D}}$, reveal how collective diffusion changes due to self-crowding and competition between subpopulations.

In what follows we extend model (2) to incorporate reactions between subpopulations. We introduce switching at the discrete level by allowing particles to spontaneously switch between the blue and red types at a constant rate that is independent of both time and space.

II.1 Particle-based model

We consider again a population of $N$ spheres that can be either tagged blue or red, all of diameter $\epsilon$ in $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, bounded and of dimensionless volume of order one. As before, we denote the position of the $i$th individual as ${\mathbf{X}}_{i}(t)$, and $\mathcal{I}_{\mathsf{b}}(t)$ and $\mathcal{I}_{\mathsf{r}}(t)$ the sets of indices corresponding to blue and red particles, respectively. However, we note that now these sets now depend on time since particles can change color. We denote by $S_{i}(t)$ the color of the $i$th individual at time $t$, where $S_{i}(t)=\mathsf{b}$ if the particle is blue or $S_{i}(t)=\mathsf{r}$ if it is red. As before particles evolve according to overdamped Langevin SDEs:

[TABLE]

where $k_{\mathsf{b}}$ and $k_{\mathsf{r}}$ are the switching rates. When considering the probabilistic description of the discrete model, the color of the $i$th particle is a random variable. We therefore consider the joint probability density function $P(\vec{x},\vec{s},t)$ $=\mathbb{P}(\{({\bf X}_{1}(t),\ldots,{\bf X}_{N}(t))=\vec{x}\}\cap\{(S_{1}(t),\ldots,S_{N}(t))=\vec{s}\})$, the probability of the $N$ particles being in the configuration $\vec{x}=({\bf x}_{1},\ldots,{\bf x}_{N})$ with the colors of the particles being described by the vector of states $\vec{s}=(s_{1},\ldots,s_{N})$ at time $t$. The configuration space (set of legal configurations) is defined as $\Omega_{\epsilon}^{N}\times\{\mathsf{b},\mathsf{r}\}^{N}$, where $\Omega_{\epsilon}^{N}=\Omega^{N}\setminus\mathcal{B}_{\epsilon}$, and $\mathcal{B}_{\epsilon}$ is the set of illegal spatial configurations,

[TABLE]

Using the Chapman–Kolmogorov equation Ross (2014) we can write the following Fokker–Planck equation as an exact description of the evolution of the density $P(\vec{x},\vec{s},t)$ in configuration space $\Omega_{\epsilon}^{N}\times\{\mathsf{b},\mathsf{r}\}^{N}$,

[TABLE]

II.2 Population-based model

Although linear, the PDE for the particle-based model (5a) is very high-dimensional for large $N$. For this reason, we seek to derive a coupled system of low dimensional PDEs for the probability distribution of the position and state of a typical particle. Accordingly we introduce

[TABLE]

The particle choice is unimportant since $P$ is invariant with respect to permutations of particle position labels. We remark that the functions $b({\bf x},t)$ and $r({\bf x},t)$ are not marginal distribution functions as in Bruna and Chapman (2012b). Here they represent the probabilities of finding a particle at position ${\bf x}$ at time $t$, and the particle being either blue or red, respectively. We denote by $N_{\mathsf{b}}(t)$ and $N_{\mathsf{r}}(t)$ the number of blue and red particles at time $t$ and note from equations (6) that $\int_{\Omega}b({\bf x},t)\mathrm{d}{\bf x}=N_{\mathsf{b}}(t)/N$ and $\int_{\Omega}r({\bf x},t)\mathrm{d}{\bf x}=N_{\mathsf{r}}(t)/N$. We shall refer to $b({\bf x},t)$ and $r({\bf x},t)$ as the blue and red density functions.

As before, we focus on the low volume fraction regime $N\epsilon^{d}\ll 1$ in which pairwise interactions dominate those involving three or more particles. Under these assumptions the problem reduces to the case when $N=2$ Bruna and Chapman (2012b). We introduce the following notation, $P_{s_{1}s_{2}}({\bf x}_{1},{\bf x}_{2},t)\equiv P({\bf x}_{1},{\bf x}_{2},s_{1},s_{2},t)$. For simplicity, we describe the derivation for ${\bf f}_{\mathsf{b}}={\bf f}_{\mathsf{r}}={\bf 0}$, but the case with drift can be obtained in a similar manner. Setting $N=2$ in equation (5a) we find

[TABLE]

where $\vec{x}=({\mathbf{x}}_{1},{\mathbf{x}}_{2})$, with no-flux boundary conditions as in (5b). The initial conditions on $P_{s_{1}s_{2}}$ follow from the initial condition on $P$. For example, $P_{\mathsf{b}\mathsf{b}}({\mathbf{x}}_{1},{\mathbf{x}}_{2},0)=P_{0}({\mathbf{x}}_{1},{\mathbf{x}}_{2},\mathsf{b},\mathsf{b})$. For $N=2$, equations (6) become

[TABLE]

Combining equations (7a), (7b) and (8a) we arrive at the following equation for $b({\bf x}_{1},t)$

[TABLE]

By using the divergence and the Reynolds transport theorems, we simplify equation (9) to give the following integro-partial differential equation

[TABLE]

II.3 Matched asymptotic expansions for the densities $P_{\mathsf{b}\mathsf{b}}$ and $P_{\mathsf{r}\mathsf{r}}$

We define two regions in configuration space, the inner region where the two particles are close to each other, $\|{\bf x}_{1}-{\bf x}_{2}\|\sim\epsilon$, and the outer region where they are far apart, $\|{\bf x}_{1}-{\bf x}_{2}\|\gg\epsilon$. To evaluate the collision integral in equation (10a) we must first calculate an asymptotic expansion for $P_{\mathsf{b}\mathsf{b}}({\bf x}_{1},{\bf x}_{2},t)$ in the inner region. In the outer region the positions of the two particles are assumed to be uncorrelated at leading order. Therefore

[TABLE]

for some function $q_{\mathsf{b}}$. In the inner region, we introduce the inner variables ${\bf x}_{1}=\hat{{\bf x}}_{1}$ and ${\bf x}_{2}=\hat{{\bf x}}_{1}+\epsilon\hat{{\bf x}}$ and define $\hat{P}_{\mathsf{b}\mathsf{b}}(\hat{{\bf x}}_{1},\hat{{\bf x}},t)=P_{\mathsf{b}\mathsf{b}}({\bf x}_{1},{\bf x}_{2},t)$. Then equation (7a) becomes

[TABLE]

We seek a regular power series expansion of the form $\hat{P}_{\mathsf{b}\mathsf{b}}=\hat{P}_{\mathsf{b}\mathsf{b}}^{(0)}+\epsilon\hat{P}_{\mathsf{b}\mathsf{b}}^{(1)}+\cdots$. Substituting in (12) and equating terms of zeroth order we obtain the following problem

[TABLE]

with solution $\hat{P}_{\mathsf{b}\mathsf{b}}^{(0)}=q_{b}^{2}(\hat{{\bf x}}_{1},t)$. Similarly, the first order problem for $\hat{P}_{\mathsf{b}\mathsf{b}}^{(1)}$ is

[TABLE]

with solution $\hat{P}_{\mathsf{b}\mathsf{b}}^{(1)}=q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)\hat{{\bf x}}\cdot\nabla_{\hat{{\bf x}}_{1}}q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)+P_{\mathsf{b}\mathsf{b}}^{\text{out},(1)}(\hat{{\bf x}}_{1},\hat{{\bf x}}_{1},t)$. We conclude that, to first order in $\epsilon$, the expansion of $P_{\mathsf{b}\mathsf{b}}({\bf x}_{1},{\bf x}_{2},t)$ in the inner region is

[TABLE]

Similarly, to first order we can write the expansion of $P_{\mathsf{r}\mathsf{r}}({\bf x}_{1},{\bf x}_{2},t)$ in the inner region as

[TABLE]

II.4 Matched asymptotic expansions for the densities $P_{\mathsf{b}\mathsf{r}}$ and $P_{\mathsf{r}\mathsf{b}}$

In order to evaluate the collision integral in (10a) we require an asymptotic expansion for $P_{\mathsf{b}\mathsf{r}}({\bf x}_{1},{\bf x}_{2},t)$ in the inner region. In the outer region, our assumption of independence enables us to write

[TABLE]

for some functions $q_{\mathsf{b}}$ and $q_{\mathsf{r}}$.

[TABLE]

As before, we seek a solution of the form $\hat{P}_{\mathsf{b}\mathsf{r}}=\hat{P}_{\mathsf{b}\mathsf{r}}^{(0)}+\epsilon\hat{P}_{\mathsf{b}\mathsf{r}}^{(1)}+\cdots$. Repeating the procedure outlined in the previous section, we find that the inner region solution for a blue and a red particle is

[TABLE]

where $q_{\mathsf{b}}$ and $q_{\mathsf{r}}$ are evaluated at $(\hat{{\bf x}}_{1},t)$. Similarly, to first order we find that the expansion of $P_{\mathsf{r}\mathsf{b}}({\bf x}_{1},{\bf x}_{2},t)$ in the inner region is

[TABLE]

II.5 Evaluating the collision integral

Since the collision integral, which we shall denote by $I$, in (10a) is defined over the boundary surface $\partial B_{\epsilon}({\bf x}_{1})$, we use the inner solutions $\hat{P}_{\mathsf{b}\mathsf{b}}$ and $\hat{P}_{\mathsf{b}\mathsf{r}}$ in order to evaluate it. We first transform $I$ into inner variables and use the boundary conditions (12b) and (16b) yielding

[TABLE]

Inserting the inner region solutions (13) and (17) into (19), we find

[TABLE]

where the constants $\omega$, $\xi_{\mathsf{b}}$ and $\eta_{\mathsf{b}}$ are given in (3). Using the normalization condition $\int_{\Omega_{\epsilon}^{2}}\left(P_{\mathsf{b}\mathsf{b}}+P_{\mathsf{b}\mathsf{r}}+P_{\mathsf{r}\mathsf{b}}+P_{\mathsf{r}\mathsf{r}}\right)\mathrm{d}{\mathbf{x}}_{1}{\mathbf{x}}_{2}=1$, we find that $q_{\mathsf{b}}({\bf x}_{1},t)=b({\bf x}_{1},t)+O(\epsilon^{2})$ and $q_{\mathsf{r}}({\bf x}_{1},t)=r({\bf x}_{1},t)+O(\epsilon^{2})$. This allows us to replace $q_{\mathsf{b}}$ and $q_{\mathsf{r}}$ in (20) by $b$ and $r$, respectively.

Substituting the asymptotic expansion for the collision integral (20) into equation (10a) we arrive at the following equation for the blue particles’ density $b({\mathbf{x}},t)$ (after replacing ${\mathbf{x}}_{1}$ by ${\mathbf{x}}$):

[TABLE]

A similar calculation on (10b) yields an equation for the density of red particles $r({\mathbf{x}},t)$.

For the case of $N$ particles there is a contribution from every potential pairwise interaction. This allows us to write down the population-level equation for the general case, where we reintroduce the drift contributions, as follows

[TABLE]

II.6 Numerical results

In this section we present numerical simulations of the discrete stochastic model (4) and the population-based model (21). The discrete model was simulated using the software package Smoldyn Andrews and Bray (2004); Andrews (2017) via the standard Euler–Maruyama method Erban and Chapman (2007).

If two particles are found to be overlapping then we incorporate volume exclusion by using Smoldyn’s ballistic collision update rule. This method assigns to each particle a post-collision velocity by comparing both particles’ current and previous positions. Then it uses conservation of momentum to compute their positions had they undergone an elastic collision.

When considering spontaneous switching of particles between subpopulations, at each timestep in the algorithm we generate a random number $\rho$ uniformly on the interval $[0,1]$. A blue particle becomes red if $\rho<k_{\mathsf{b}}\Delta t$. Similarly, a red particle becomes blue if $\rho<k_{\mathsf{r}}\Delta t$. In order for this implementation to be accurate we assume $\Delta t$ to be small enough such that $k_{\mathsf{b}}\Delta t\ll 1$ and $k_{\mathsf{r}}\Delta t\ll 1$, and so are good approximations to the probability of a switch occurring in a time interval of length $\Delta t$. We also require that the timestep is small enough that each particle moves, on average, a distance much less than the particle’s diameter. This requirement ensures that most collisions are detected.

We solved the PDEs numerically using the method of lines. In this technique only the spatial variables are discretized and time is viewed as a continuous variable, resulting in a system of Ordinary Differential Equations (ODEs) that are solved using MATLAB’s inbuilt ode15s solver (which performs the discretization of the time variable adaptively).

Initially we have $350$ blue particles and $50$ red particles in the domain $\Omega=[-1/2,1/2]^{2}$. We select an initial Gaussian density of particles with zero mean and a standard deviation of $0.09$ for both the blue and red subpopulations. We set $D_{\mathsf{b}}=0.5$ and $D_{\mathsf{r}}=1$ for the diffusion coefficients of the blue and red populations respectively, as well as switching parameters $k_{\mathsf{b}}=k_{\mathsf{r}}=10$. We ignore drift terms, setting ${\bf f}_{\mathsf{b}}={\bf f}_{\mathsf{r}}={\bf 0}$. The data from the stochastic simulations of the SDEs (4a) and (4b) were generated using $10^{4}$ realizations, which resulted in $4\times 10^{6}$ individual trajectories, where a timestep of $\Delta t=10^{-5}$ was used. The stochastic simulations used to generate the results presented in this section were performed on an AMD FX(tm)-4350 CPU, with an advertised processor speed of $4.2$ GHz. The total CPU time needed to generate the trajectories seen in Figure 1 is approximately $120$ hours, whereas the PDEs (21) are solved numerically within $13$ seconds.

Figure 1 shows the densities $b({\mathbf{x}},t)$ and $r({\mathbf{x}},t)$ of the blue and red particles respectively at time $t=0.05$. Panels (I) and (II) respectively show the PDE solution and the histogram produced from the simulated trajectories for point particles ($\epsilon=0$), whereas panels (III) and (IV) respectively show the same for finite size particles of diameter $\epsilon=0.01$. Figure 2(a) shows a one dimensional slice of the histograms in Figure 1 for a clearer comparison of the PDEs and the stochastic data. We observe good agreement between the particle-based model and the derived population-level equations. The relative $L_{2}$-norm errors in the densities (between the solution of the PDE model and the histograms of the stochastic particle-based model) are, from top to bottom in Figure 2(a), $2.5\%,1.7\%,2.6\%$ and $2.7\%$.

For further quantitative comparison we consider how the population numbers, $N_{\mathsf{b}}(t)$ and $N_{\mathsf{r}}(t)$, evolve over time. Recall that $\int_{\Omega}b({\bf x},t)\mathrm{d}{\bf x}=N_{\mathsf{b}}(t)/N$ and $\int_{\Omega}r({\bf x},t)\mathrm{d}{\bf x}=N_{\mathsf{r}}(t)/N$. We can then integrate equations (21a) over the domain $\Omega$ to obtain, using zero-flux boundary conditions, the following ODE

[TABLE]

Figure 2(b) shows the solution of (22) and compares them with results from the stochastic simulations. We only plot results for finite size particles but note that the curves for point particles are identical. That is, the finite size particles do not effect the evolution of the population number when the switching rates are constant on $\Omega$, and so the ODEs are exactly the same as if we had considered the non-spatial model originally.

II.7 Application to a hybrid model of chemotaxis

We now consider a three species model of chemotaxis as an application of our proposed framework. As before we consider two populations of finite-size particles, blue and red, which we shall refer to as cells. We introduce a third species which shall act as a chemoattractant that is produced by the red cells and acts as a signaling chemical towards which the blue cells migrate. The size of chemoattractant molecules is negligible in comparison to the size of the cells which produce them. As such we consider the chemoattractant as a continuum field described by a reaction-diffusion PDE. This yields a hybrid modeling framework, whereby some species are described as a continuum and some species as a set of particles Franz and Erban (2013). Hybrid chemotaxis models have been developed in the literature Dallon and Othmer (1997); McLennan et al. (2012). In Dallon and Othmer (1997), a population of identical cells that only interact via the chemical signal (that is, cells are represented as point particles) is considered. The model in McLennan et al. (2012) describes cell migration during embryonic development and includes two types of cells (“leaders” and “followers”). Cells undergo a position-jump process (with jumps of a fixed distance) and volume exclusion between cells is included as a rejection mechanism (cells are modeled as hard discs and a jump is aborted if a cell attempts to move to a position where it overlaps with another cell).

In this section, we use our model for two reacting species of hard-sphere particles in the context of hybrid modeling of chemotaxis. In contrast to McLennan et al. (2012), our model is for Brownian particles and includes a more accurate description of excluded-volume interactions (not biasing the motion of cells by aborting moves that would lead to overlaps Strating (1999)). The red cells secrete chemoattractant at a constant rate $\gamma$, while the blue cells undergo biased motion up spatial gradients of the chemoattractant. The chemoattractant diffuses throughout the domain without interacting with the cells and degrades at a constant rate $\mu$. An additional mechanism for the loss of chemoattractant due to consumption by blue cells could be considered here. However, under the assumption that the dominant mechanism for the loss of chemoattractant is natural decay Grima (2005), the additional consumption term is neglected. We assume cells move in a two-dimensional domain $\Omega$ and that the chemical is allowed anywhere in the domain (to represent cells in contact with the substrate and the chemical diffusing in the space above). We denote the position and the color of the $i$th cell at time $t$ as ${\bf X}_{i}(t)$ and $S_{i}(t)\in\{\mathsf{b},\mathsf{r}\}$ respectively. We denote by $C({\bf x},t|\vec{x},\vec{s})$ the concentration of the chemoattractant at position ${\bf x}\in\Omega$ at time $t$ given a configuration state $\vec{x}=\{{\bf X}_{1}(t),\ldots,{\bf X}_{N}(t)\}$ and $\vec{s}=\{S_{1}(t),\ldots,S_{N}(t)\}$. The microscopic model for the blue cells ($S_{i}(t)=\mathsf{b}$) is

[TABLE]

Finally, cells change color according to the reactions

[TABLE]

Blue cells switch color at a rate proportional to the concentration of chemoattractant at their current position $X_{i}(t)$, and red cells switch color at a constant rate $k_{\mathsf{r}}$ as before. Similar to Section II.1, we can formulate the microscopic model in terms of a high-dimensional Fokker–Planck equation and introduce the population-level equations $b({\bf x},t)$ and $r({\bf x},t)$ as in (6). Additionally introducing the population-level equation for the chemoattractant $c({\bf x},t)=\int_{\Omega_{\epsilon}^{N}}\sum_{\vec{s}\in\mathcal{C}^{N}}C({\bf x},t|\vec{x},\vec{s})\mathrm{d}\vec{x}$ allows us to write down the following population-level reaction-diffusion equations

[TABLE]

As before, this model includes the excluded-volume effects up to $O(\epsilon^{2})$ in the equations for $b$ and $r$. Regarding the chemical, the error in the reduction from the kernel $K_{h}$ in the equation for $C$ to the source term $\gamma Nr$ in the equation for $c$ is $O(h^{2})$ (this can be seen using, for example, Laplace’s method). Hence our choice of $h$ ensures that this error is smaller than $O(\epsilon^{2})$.

We conclude this Subsection by performing numerical simulations of the hybrid model for chemotaxis (23) and the reaction-diffusion PDE system (24). Since a key feature of this model is the chemotactic drift, we must adapt our stochastic simulations to account for the drift in the SDEs of blue particles (23a). Since Smoldyn Andrews (2017) does not presently allow for biased Brownian motion, we perform the particle-based simulations in this Subsection using the alternative library Aboria Robinson and Bruna (2017); Robinson (2018). The SDEs (23a) and (23b) are again integrated using the Euler–Maruyama method with a constant timestep $\Delta t$. The same timestep is used to integrate the reaction-diffusion PDE (23c), using the Euler method in time and a second-order finite-difference scheme in space, and to simulate the reactions (23d). For details on the implementation of the hybrid model with Aboria we refer the readers to Bruna et al. (2018).

Figure 3 shows a comparison between the hybrid model (23) and continuum model (24) of chemotaxis with and without excluded-volume interactions. We start with a population of $N=800$ cells of diameter $\epsilon=0.01$ in the domain $\Omega=[-1/2,1/2]^{2}$. Half of the particles are initially blue and uniformly distributed, and the other half are red and uniformly distributed in a vertical band centered at $x=0$ of width 3/31. The initial condition for the chemical is $c_{0}(x)=1$. We perform 5000 realizations of the stochastic simulation so that, as before, we have $4\times 10^{6}$ trajectories. Figures 3(a-c) show the densities $b$, $r$, and the concentration $c$ averaged over the vertical direction at four different times $t$. As time progresses, the red particles produce chemoattractant and this leads to a higher concentration of chemical in the centre of the domain (see Figure 3(c)). This gradient in concentration induces a drift to the blue particles, that want to move towards the middle of the domain. In the case of point particles, this gives rise to a hill-shaped profile of the blue density centered at the origin. However, in the case of finite-size particles, the interplay between the gradient of concentration $\nabla_{\mathbf{x}}c$ and the gradient of reds $\nabla_{\mathbf{x}}r$ results in a volcano shaped profile of the blue density (see Figure 3(a)). There is a net production of red particles to the detriment of the blue particles and a net decay of chemical $c$ (Figure 3(d)). We note that, as for the example in Figure 2(b), the time evolution of the total number of blue and red particles is not affected by the excluded-volume interactions. We find good agreement between the hybrid and continuum models. The relative $L_{2}$-norm errors in the densities in Figures 3(a-c) (using the solutions at the histogram grid points and $t=0,0.01,0.02,\dots,0.05$) are, for point particles ($\epsilon=0$), $0.39\%,0.85\%$, and $0.038\%$, respectively. The errors in the case of $\epsilon=0.01$ are $0.80\%,4.2\%$, and $0.15\%$, respectively. We note that these errors are very small compared to errors we would have committed by ignoring the excluded-volume effects: the errors between the PDE model for $\epsilon=0$ and the stochastic model for $\epsilon=0.01$ are $12\%,13\%$ and $0.40\%$ for the blue particles, red particles, and chemical concentration, respectively.

III Switching Upon Collision: Bimolecular reactions

In this section we assume that particles switch between subpopulations as a result of a pairwise collision or interaction and consider bimolecular reactions of the form $B+R\rightarrow 2R$. We first describe how this effect changes our particle-level description and then derive the corresponding population-level equations.

III.1 Particle-based model

Our discrete model has the same Langevin dynamics as before (4). When two blue (red) particles are in contact, $\|{\bf X}_{i}-{\bf X}_{j}\|=\epsilon$ and $i,j\in\mathcal{I}_{\mathsf{b}}(t)$ or $\mathcal{I}_{\mathsf{r}}(t)$, we impose standard reflective boundary conditions. However, when a blue and a red particle are in contact, $i\in\mathcal{I}_{\mathsf{b}}(t)$ and $j\in\mathcal{I}_{\mathsf{r}}(t)$, the blue particle becomes red with a certain probability, resulting in $S_{i}(t)=\mathsf{r}$ and $N_{\mathsf{b}}(t)=N_{\mathsf{b}}(t)-1$ after the collision.

As before the joint probability density function $P(\vec{x},\vec{s},t)$ obeys the linear Fokker–Planck equation

[TABLE]

for $\vec{x}\in\Omega_{\epsilon}^{N}$. The boundary conditions now require more careful treatment. We define a flux for each configuration of particle states

[TABLE]

In the case of spontaneous switching on the internal boundaries $\partial\Omega_{\epsilon}^{N}\setminus\partial\Omega^{N}$ we had that $\mathcal{J}_{\vec{s}}=0$ for all $\vec{s}\in\mathcal{C}^{N}$. However there are now non-zero fluxes across these internal boundaries in the configuration space, where the reactivity parameter for these boundaries is $\lambda$. To make this explicit, we consider the boundary between a pair of particles, $\|{\bf x}_{i}-{\bf x}_{j}\|=\epsilon$ for $i\neq j$. If ($s_{i},s_{j}$)=($\mathsf{b},\mathsf{b}$) then $\mathcal{J}_{\vec{s}}=0$ on this boundary, as there is no reaction between two blue particles. If ($s_{i},s_{j}$)=($\mathsf{b},\mathsf{r}$) or ($\mathsf{r},\mathsf{b}$), then we have $\mathcal{J}_{\vec{s}}=\lambda P_{\vec{s}}$ on $\|{\bf x}_{i}-{\bf x}_{j}\|=\epsilon$. This is because there is a flux out of the configuration state $\vec{s}$ due to the bimolecular reaction. Finally if ($s_{i},s_{j}$)=($\mathsf{r},\mathsf{r}$) no reaction occurs for two red particles; however there is a flux into the configuration state $\vec{s}$ on the boundary $\|{\bf x}_{i}-{\bf x}_{j}\|=\epsilon$ from the blue and red reactions. Therefore on this boundary $\mathcal{J}_{\vec{s}}=-\lambda\left(P_{\vec{c}_{i}}+P_{\vec{c}_{j}}\right)$, where $\vec{c}_{i}=(s_{1},\ldots,s_{i-1},\mathsf{b},s_{i+1},\ldots,s_{N})$ and $\vec{c}_{j}=(s_{1},\ldots,s_{j-1},\mathsf{b},s_{j+1},\ldots,s_{N})$.

III.2 Population-based model

To reduce the particle-level description to the population-level we consider as before the density functions (6). We restrict attention to the case for $N=2$ and rewrite the problem in terms of $P_{s_{1}s_{2}}({\bf x}_{1},{\bf x}_{2},t)=P({\bf x}_{1},{\bf x}_{2},s_{1},s_{2},t)$ as before. For simplicity, we again set ${\bf f}_{\mathsf{b}}={\bf f}_{\mathsf{r}}={\bf 0}$. For $N=2$ we only have one inner boundary, $\|{\bf x}_{1}-{\bf x}_{2}\|=\epsilon$, with the following boundary conditions

[TABLE]

Then, from (25) the equations for the blue and red densities with the collision integrals can be written as

[TABLE]

We next compute the asymptotic expansions for $P_{\mathsf{b}\mathsf{b}}$, $P_{\mathsf{b}\mathsf{r}}$, $P_{\mathsf{r}\mathsf{b}}$ and $P_{\mathsf{r}\mathsf{r}}$ in the inner region in order to evaluate the integrals in (28).

III.3 Matched asymptotic expansions and collision integral

First we consider $P_{\mathsf{b}\mathsf{b}}({\bf x}_{1},{\bf x}_{2},t)$. As the boundary condition on $\|{\bf x}_{1}-{\bf x}_{2}\|=\epsilon$ is $\mathcal{J}_{\mathsf{b}\mathsf{b}}=0$, we can use the results of Section II.3 to write down the inner solution as

[TABLE]

for some function $q_{\mathsf{b}}$.

We construct the inner solution to $P_{\mathsf{b}\mathsf{r}}({\bf x}_{1},{\bf x}_{2},t)$ by matching the solutions in the outer and inner regions. Introducing the inner variables ${\bf x}_{1}=\hat{{\bf x}}_{1}$ and ${\bf x}_{2}=\hat{{\bf x}}_{1}+\epsilon\hat{{\bf x}}$ as before, we transform equations (25) and (27) to get the following problem

[TABLE]

on $\|{\bf x}\|=1$, and matching condition given by (16c) as $\hat{\mathbf{x}}\to\infty$. We seek a regular power series expansion of the form $\hat{P}_{\mathsf{b}\mathsf{r}}=\hat{P}_{\mathsf{b}\mathsf{r}}^{(0)}+\epsilon\hat{P}_{\mathsf{b}\mathsf{r}}^{(1)}+\cdots$, and as in Section II the zeroth order problem has the trivial solution $\hat{P}_{\mathsf{b}\mathsf{r}}^{(0)}=q_{\mathsf{b}}(\hat{{\bf x}}_{1},t)q_{\mathsf{r}}(\hat{{\bf x}}_{1},t)$. From (30), the problem for $\hat{P}_{\mathsf{b}\mathsf{r}}^{(1)}$ is

[TABLE]

In two dimensions ($d=2$), problem (31) does not have a solution satisfying both the inner and outer boundary conditions, indicating that a more intricate method is required. For the case $d=3$ the solution is

[TABLE]

Combining the results above we find that

[TABLE]

We denote by $I$ the collision integral in equation (28a),

[TABLE]

Using the inner solutions (29) and (32) we can now evaluate $I$. Expressing it in terms of inner variables and using the boundary condition (30b), we obtain

[TABLE]

where the constants $\omega$, $\xi_{\mathsf{b}}$ and $\eta_{\mathsf{b}}$ are defined in (3). We use again the normalization condition to replace $q_{\mathsf{b}}$ and $q_{\mathsf{r}}$ by $b$ and $r$, respectively. Combining these results we deduce that the population-level equation for a blue particles $b({\mathbf{x}},t)$ satisfies

[TABLE]

From (28b) a similar equation can be derived for the red density. For the general case of $N$ particles, reintroducing drift terms, we have the following system of nonlinear reaction-diffusion-advection equations

[TABLE]

for ${\bf x}\in\Omega$, where the diffusion matrix, $\mathcal{D}$, and the drift matrix, ${\mathcal{F}}$, are defined by equations (21b) and (21c). The reaction terms $\mathcal{R}$ are given by

[TABLE]

We interpret the reaction terms as follows: $(N-1)4\pi\epsilon^{2}$ is the total surface area of the $(N-1)$ spheres with which the sphere under consideration could react. The factor $\left(1-\epsilon\lambda/(D_{\mathsf{b}}+D_{\mathsf{r}})\right)$ represents a correction to the source and sink terms, and indicates that only a fraction of the available surface area is used. If the diffusion coefficients are large relative to the reactivity parameter $\lambda$ then this correction term approaches one, which can be interpreted as the case when the particle is moving fast enough that the entire surface area is available to it for collisions.

Our result (36) can be related to the $\lambda$ - $\bar{\rho}$ model introduced by Erban and Chapman Erban and Chapman (2009) as follows. In their model, a bimolecular reaction can only take place when the distance between two particles is less than $\bar{\rho}$, and it occurs at a rate $\lambda$. They distinguish between two cases depending on whether the relative mean-square displacement between the particles in one timestep, $\sqrt{2(D_{\mathsf{b}}+D_{\mathsf{r}})\Delta t}$, is small or large compared to the reaction radius $\bar{\rho}$. We consider the former case, since we already perform simulations in this regime so that we do not miss collisions between particles (as discussed in Section II.6). Taking the limit $\lambda$ small, the reaction rate was found to be proportional to the volume of the reactive region (see Eq. (31) in Erban and Chapman (2009)). By contrast, the leading order of our reaction rate is proportional to the surface area of the collision surface (see (36)). This discrepancy arises because we have incorporated the excluded-volume effect directly into our microscopic model and, further, particles are not allowed to get closer to each other than the reaction radius, that is, $\bar{\rho}\equiv\epsilon$ in our case.

III.4 Numerical Results

In this section we show how the inclusion of bimolecular reactions affects the time dependent behavior of our model. The stochastic simulations are implemented as described previously in Section II.6 with the exception of the reaction terms. When considering particle switching due to the collisions we proceed as follows. When a red and blue particle overlap, the reaction is implemented by generating a random number $\rho$ uniformly on the interval $[0,1]$. If $\rho<\mathcal{P}\sqrt{\Delta t}$, where $\mathcal{P}\sqrt{\Delta t}$ is the reaction probability, we switch the color of the blue particle to red and then update the positions due to the collision. We note that the relationship between the reactivity parameter, $\lambda$, and the reaction probability parameter, $\mathcal{P}$, has been shown to be Erban and Chapman (2007)

[TABLE]

In order to increase the accuracy of this method we may also want to account for the possibility, that during the interval $[t,t+\Delta t]$, the red and blue particles may have overlapped and then moved apart, and a reaction may have occurred that otherwise would be missed. This correction, first pointed out by Andrews and Bray Andrews and Bray (2004), is not included in our implementation of the bimolecular reactions. Instead, the timestep $\Delta t$ is taken to be small enough that the probability of particles overlapping and separating during $\Delta t$ is negligible Andrews and Bray (2004).

We consider the evolution of the population numbers by integrating equations (35) over the domain $\Omega$ to obtain the following integro-differential equation for the number of blue particles $N_{\mathsf{b}}(t)$

[TABLE]

and $N_{\mathsf{r}}(t)=N-N_{\mathsf{b}}(t)$. Since the leading-order reaction term in equation (38) is at order $\epsilon^{2}$, it is natural to ask if retaining terms to $O(\epsilon^{2})$ provides sufficient detail to accurately predict the evolution of the population numbers. In Figure 4(a) we present the solutions to equation (38) correct to $O(\epsilon^{2})$ and $O(\epsilon^{3})$. We note that the $O(\epsilon^{3})$ terms are needed to accurately model the evolution of the population number. In particular, the relative $L_{2}$-norm errors of the PDE model with only the $O(\epsilon^{2})$ reaction terms are $5.9\%$ and $4.0\%$ for blue and red particles respectively. The relative errors when adding the $O(\epsilon^{3})$ reactions terms in the PDE are reduced to $0.9\%$ and $0.6\%$, respectively.

We also consider how the total volume fraction occupied by spheres affects the evolution of the population of particles. To investigate the effect of volume fraction we plotted the evolution of the blue population for a variety of different diameters $\epsilon$ in Figure 4(b), keeping $N$ fixed. We see that the time to extinction decreases as volume fraction increases. In particular, when only $6.55\%$ of the domain is filled with particles we observe rapid extinction of the blue particles population.

IV Discussion

In this paper we have studied the interplay between volume exclusion and reactions in a two-species system of Brownian hard-sphere particles. Reactions between the two subpopulations of particles, which we termed as either “blue” or “red” particles, are either spontaneous or as a result of hard-core collisions.

In Section II we considered the case of spontaneous switching. We formulated the discrete model, a system of overdamped Langevin SDEs, to allow particles to switch their color at a spatially dependent rate. The resulting system of nonlinear reaction-diffusion-advection equations (21) is similar to that for multiple populations without switching Bruna and Chapman (2012b) with additional linear source and sink terms. We then presented numerical solutions of the PDEs and showed that they are in good agreement with stochastic simulations of the discrete model. We also showed how our modeling framework can be adapted to a well-studied problem in mathematical biology, namely that of cell chemotaxis. In this case, the spontaneous reaction from blue to red particles is mediated by the chemical. By introducing a third species of point particles representing the chemoattractant, we showed how to coarse-grain a hybrid model whereby cells are modeled as individual particles of a finite-size and the chemical as a continuum already, to a continuum model for all three species. The hybrid model is similar to that used by McLennan et al. (2012), but here we were able to show how the hybrid model upscales to a continuum model. In particular, numerical simulations of the hybrid and continuum models showed that excluded-volume interactions in combination with reactions and chemotaxis can have a significant impact on the distribution of cells in the domain.

In Section III we considered bimolecular reactions due to collisions between the hard-sphere particles. In particular, we considered the reaction $B+R\rightarrow 2R$, so that a blue particle may change to red, with a prescribed probability, after a collision with a red particle. This bimolecular reaction comes into the model as a nontrivial boundary condition in the high-dimensional Fokker–Planck PDE and, as a result, it changes the problem to solve via matched asymptotic expansions substantially. We obtained a solution for the three-dimensional case. The case of two spatial dimensions is harder, because a random walker will return to its starting point with probability one in two dimensions, but not in three Durrett (2013). We leave the two-dimensional case, which could be solved by introducing an intermediate region between the inner and outer regions in the method of matched asymptotic expansions, for future work. The cross-diffusion model with bimolecular reactions between subpopulations is the main contribution of this paper. To our knowledge, this is the first continuum PDE model for (off-lattice) Brownian particles that combines diffusion and reactions occurring due to contact between particles, instead of introducing an artificial reaction radius Erban and Chapman (2009).

The method presented in this paper provides a systematic way of linking the stochastic particle-level description and the population-level PDE description. As such it is ideally suited for experimental validation. Depending on the type of experimental data available for a given biological system, either the stochastic model or the PDE model might be more readily parameterized. Then one could exploit the link between the two levels of description to validate the method and learn more about the experimental system. This study is beyond the scope of the present work.

An interesting extension of this framework would be to incorporate cell proliferation and death. The main challenges when including reactions of this type into the modeling framework is that the total number of cells in the system, $N(t)$, becomes a random variable and as such the dimensionality of the configuration space can fluctuate. A further consideration is that of the volume fraction occupied in the domain: unless parent cells divide to produce offspring (keeping the total volume constant), the volume occupied in the domain will increase.

Acknowledgements.

The authors thank Dr. Steve Andrews and Dr. Martin Robinson for helpful discussions on the implementation of the stochastic algorithms using the software libraries Smoldyn Andrews (2017) and Aboria Robinson (2018), respectively. DBW would like to thank the EPSRC (Grant number EP/G03706X/1) for funding through a studentship at the Systems Biology programme of The University of Oxford’s Doctoral Training Centre.