1 Introduction
Biology has recently experienced a revolution through the development of new, quantitative, measurement techniques. As an example, superresolution microscopy allows biomolecules to be localised, counted and distinguished
(Betzig et al, 2006). Traditional modelling approaches using partial differential equations (PDEs) for populationlevel variables (concentrations, densities) are unable to capture and explain some of the particlelevel (for example, molecule or celllevel) features that we are now obtaining from experiments. This is why a more detailed modelling approach, namely particlebased or agentbased models that describe the behaviour of each organism, is now increasingly used by the mathematical biology community. Computer simulations are ideal for studying the dynamical mechanisms arising from the interplay of these particles, filling in the details that cannot be resolved experimentally, and testing/generating biological hypothesis. These simulations comprise algorithms for particlebased stochastic reactiondiffusion processes that track individual particles and are computationally expensive.The computational requirements of particlebased models, such as the efficient calculation of interactions between particles, are challenging to implement in a way that scales well with the number of particles, uniform and nonuniform particle distributions, different spatial dimensions and periodicity. The majority of existing software is typically designed to suit a particular application, and therefore tends not to be sufficiently generic that can be used to implement these computationally challenging routines, and so they are reimplemented again and again in each software package.
The primary class of software for particlebased methods is geared towards molecular dynamics simulations, such as the hugely popular GROMACS (Abraham et al, 2015), or OpenMM (Eastman et al, 2012) packages. These generally include sophisticated neighbour searching and evaluation of longrange forces, but do not include the possibility of reactions between particles, a vital component of many biological models. On the other hand, the biochemical simulator package Smoldyn (Andrews and Bray, 2004) implements both unimolecular and bimolecular reactions using the Smoluchowski method, but has limited capabilities to implement interactions between particles and transport mechanisms other than unbiased Brownian motion.
Aboria is a C++ library for the numerical implementation of particlebased models (Robinson and Bruna, 2017). Aboria aims to provide a general purpose library that allows the user complete control to define whatever interactions they wish on a given particle set, while implementing efficiently the difficult algorithmic aspects of particlebased methods, such as neighbourhood searches and the calculation of longrange forces. It is agnostic to any particular numerical method, for example it does not contain any particular molecular dynamics or Smoluchowski dynamics algorithms, but instead provides the user with computational tools that allow them to implement these algorithms in a customised fashion. It is thus especially suitable for the implementation of novel particlebased methods, or hybrid models that couple particlebased models with continuum PDE models. In particular, Aboria gives the user:

A container class to store a particle set with a position and unique id for each particle, as well as any number of userdefined variables associated to each particle with arbitrary types (for example, to store a velocity, a force, or the particle size).

The ability to embed each particle set within an dimensional hypercube with arbitrary periodicity. The underlying data structure can be a cell list, kdtree or hyper octtree. The data structure can be chosen to suit the particular application.

Flexible neighbourhood queries that return iterators, and can use any integer norm distance measure (for example, the Euclidean distance is given by , the Chebyshev distance is given by ).

The ability to calculate longrange forces between particles efficiently (that is, in ) using the blackbox fast multipole method (Fong and Darve, 2009). This method can be used for any longrange force that is wellapproximated by a loworder polynomial for large distances between particles (the computational cost of the model scales with the degree of the polynomial).

An easy to use embedded Domain Specific Language (eDSL) expressing a wide variety of particleparticle interactions. For specialised cases which cannot be easily expressed by this symbolic eDSL, Aboria also provides a lowlevel interface based on the C++ Standard Template Library (Stroustrup, 2013).
The aim of this paper is to demonstrate how Aboria can be used to implement many features of particlebased models common in mathematical biology. These include heterogeneous particle populations, interactions between particles, volume exclusion, biased transport, reactions, and proliferation. We showcase the many features of the Aboria library through a wellknown model in mathematical biology, namely cell chemotaxis. In Sect. 2 we detail the mathematical model for chemotaxis. Then in Sect. 3 we show a stepbystep implementation of the model with Aboria. Finally, we show how to output and analyse the data using Python in Sect. 4.
2 Model for chemotaxis
As our guiding example, we consider chemotactic cell migration with volume exclusion. In particular, we study a population of cells divided into two subpopulations (types and ) and a diffusing attractive chemical substance that is produced by cells of type and consumed by cells of type . Cells move around in a twodimensional domain due to Brownian motion, biased by gradients in the concentration of the chemical in the case of type cells. Cells can also undergo reactions that change their type from chemical producers to chemical consumers. Volume exclusion can be modelled by either point particles with a shortrange repulsive interaction potential (so that cells are allowed to deform if in close contact) or hardsphere particles that are not allowed to overlap. To make the example as broad as possible, we take a hybrid modelling approach whereby the cells are modelled discretely using particlebased models for Brownian motion, while the chemical is represented by its concentration using a reactiondiffusion PDE. Hybrid models of chemotaxis have been considered in Guo et al (2008); Dallon and Othmer (1997); Franz and Erban (2013); Newman and Grima (2004); McLennan et al (2012).
Let denote the position of the th particle in . Let denote the set of cells of type and denote the set of cells of type at time . For each particle, the motion through space is described by a stochastic differential equation (SDE). For ,
(1a) 
where is the diffusion coefficient of cells of type , and denotes the gradient with respect to . For
(1b) 
Here is the diffusion coefficient of cells of type , denotes the chemotactic sensitivity of the cells (taken to be constant) and is the concentration of chemical at position and time . The interaction potential may be a soft potential incorporating effects such as size exclusion by cells and cellcell adhesion. Typical examples are a Morse potential (D’Orsogna et al, 2006; Middleton et al, 2014), an exponential potential (Bruna et al, 2017), or a LennardJones potential (Jeon et al, 2010). Alternatively, the interaction between cells may be modelled as a singular hardsphere potential so that cells are not allowed to overlap each other, for example assuming cells have diameter and taking for , otherwise would impose that for all (Bruna and Chapman, 2012). Also, cells can change type according to the following reactions:
(1c) 
Cells of type change to type with constant rate , whereas cells of type change type with rate , where is the position of the type cell (that is, the rate may depend on the chemical concentration at the location of the cell).
Finally, cells of type secrete chemoattractant at a constant rate , while cells of type consume it at a constant rate . The chemical diffuses with a diffusion constant and degrades with rate . The chemical concentration evolves according to the PDE
(1d) 
where and denote the random measures for the density of cells of type and , respectively,
(1e) 
The SDEs (1a) and (1b) and the PDE (1d) are complemented with suitable initial conditions and either periodic or noflux boundary conditions).
The numerical implementation of model (1) is challenging for several reasons. In the remainder of this section we go through each of the problems one faces. In addition to the points below, it is worth keeping in mind that generally we are interested in statistical averages of the simulations (in order to compare them with, for example, experiments or continuum PDE models). As a result, it is crucial that the simulation method is implemented efficiently, exploiting parallelisation and other algorithms to speed up the simulation.
2.1 Timestepping for diffusion and reactions
The standard way to numerically integrate the SDEs (1a) and (1b) is to use a fixed timestep and a Euler–Maruyama discretisation. For (1a), this reads
(2) 
where
is a twodimensional normally distributed random variable with zero mean and unit variance. Reactions (
1c) can also be simulated using a fixed timestep. For example, if is the number of type cells at time , then the first reaction occurs during if , whereis a uniformly distributed random number,
.The downside of the fixed timestepping approach to simulating reactions is that: (i) the timestep must be chosen to ensure that , which imposes a restriction on the size of ; (ii) choosing such a small means that in most timesteps no reactions will take place. Hence, lots of random numbers need to be generated before the reaction takes place (Erban et al, 2007); (iii) there is an exact and more efficient simulation algorithm, the Gillespie algorithm (Gillespie, 1977). The Gillespie algorithm computes the time that the next reaction will occur as , where again.
In latticebased models for reactiondiffusion processes, space is discretised into a regular lattice and diffusion is represented as jumps between neighbouring lattices, and can be treated a reaction events. This implies that both reactions and diffusion can be implemented using the same framework, the Gillespie algorithm being the obvious choice. In contrast, offlattice Brownian motion models such as (1a) and (1b) do not fit this framework and are naturally implemented with a fixed timestep approach (2). As a result, the simulation algorithm is either of fixed timestep for both processes, where is small enough to resolve diffusion, reactions, and interactions well (see Subsec. 2.2), or fixed timestep for the cell position updates and variable timestep for the cell number updates.
2.2 Cellcell interactions
Pairwise interactions between particles will generally lead to an loop to compute the interaction terms in (1a) and (1b) at every timestep . If the interaction potential is short ranged, it is convenient to use neighbourhood searches as we only need to evaluate the distances and forces between particles that are close enough. This reduces the computational cost to , where is the typical number of particles in the neighbourhood of one particle.
For longrange forces, neighbourhood searches do not help as every particle interacts significantly with every other particle. However, in many commonly used interactions forces (for example, electrostatics, gravitational) the interactions between wellspaced clusters of particles can be efficiently approximated by means of the fast multipole method (FMM) (Greengard and Rokhlin, 1987), which also leads to a total computational cost of . Aboria implements a version of the blackbox FMM (Fong and Darve, 2009)
, which uses Chebyshev interpolation to approximate the interaction of wellseparated clusters. Since the present chemotaxis model uses shortrange interactions (to represent cell volume exclusion), we do not discuss the FMM further in this paper. For more information of Aboria’s FMM capabilities, the reader is referred to the
documentation.Interactions between particles also require a careful choice of so that they are well resolved. If is too large, interactions between Brownian steps might be missed. A good rule of thumb is that the mean relative displacement (ignoring any drift terms) between two particles with diffusion coefficients and respectively, , should be less than the range of the interaction potential (equal to the sum of the particles’ radii in the case of a hardsphere interaction). It is not uncommon for shortrange potentials to be singular or very steep at the origin. In these cases, the scheme (2) may not resolve well the original SDE (1a) unless the timestep is prohibitively small (so that the drift term in (2) does not send particles very far apart, possibly missing other interactions on the way). An alternative to the explicit Euler–Maruyama scheme (2) is to use an implicit scheme with better convergence properties, or the socalled tamed Euler scheme (Hutzenthaler et al, 2012), which modifies the drift term of (1a), , where , so that it is uniformly bounded by one:
(3) 
This scheme coincides with the Euler–Maruyama scheme (2) up to order and it is just as simple to implement, but it has the advantage of allowing larger simulation timesteps for (1a) with repulsive potentials singular at the origin.
In the case of a hardsphere interaction potential, there are several options to implement the collisions between particles (see Bruna, 2012, p. 33). One option is to update particle positions according to (2) and correct any overlap at time by moving particles apart in the direction along the line joining the two particle centres. Namely, if the distance between two particles is , where is the diameter of the th particle, then the particles are moved apart a distance , where is the distance that particles have penetrated each other illegally. Note that we use twice this distance to account for the fact that particles would have collided and moved apart by (moving them apart only by would make them be exactly in contact). The way the total update distance is distributed among particles depends on the mean travelled distance of each particle and their diffusion coefficients. If the particles are of the same type, then the distance is shared equally among them. If instead one particle is immobile (suppose it is a fixed obstacle), then the total displacement is imposed on the other particle. In general, if the particles have diffusion coefficients and , the th particle takes of the displacement, and the rest goes to particle .
Another method to implement hardsphere collisions is known as the elastic collision method (Scala et al, 2007). It consists of an eventdriven algorithm between Brownian timesteps of length , whereby each Brownian particle is attributed a “velocity” and collisions between all particles in the interval are predicted and treated using a standard eventdriven method for ballistic dynamics. On the one hand, this method predicts rather than corrects collisions, and it is therefore more accurate than the first one. This implies that one may take larger steps . On the other hand, the eventdriven method is computationally more intensive and is more complex to implement.
2.3 Spatial matching between discrete and continuous variables
The hybrid model (1) combines two modelling frameworks: a particlebased approach for the cells, equations (1a) and (1b), and a continuum approach for the chemical concentration, (1d). This means that we have to consider each part of the model separately and establish a way to match them in space.
One approach is that taken by Newman and Grima (2004), where the chemical concentration is found by formally integrating equation (1d) along the cell paths. Assuming , the result is
(4) 
where and are given in (1e) and is the Green’s function for the chemical diffusion equation in ,
(5) 
Since we have an explicit expression for the chemical concentration in the whole space, the chemical gradient in (1b) and the chemical concentration in the reaction rate in (1c) can be evaluated exactly. In particular, the gradient in (1b) can be written as
(6) 
Therefore, this approach requires the history of the cell positions only, and no explicit evaluation of the chemical concentration (unless is required as a simulation output, in which case one uses (4)). One major drawback of this approach is that it only works when the domain is the whole space, and therefore it may not be applicable in many cases.
The alternative approach is to discretise and integrate the chemical concentration on a grid, for example using a finitedifferences method. This approach works for bounded domains, but it has the disadvantage that it requires spatial matching between the discrete and continuum variables. Specifically, the particles can be positioned at an arbitrary point inside the domain, while the chemical concentration is only calculated at grid points . Then we have a twoway matching to do: interpolate the concentration at the offgrid particle positions to simulate (1b) and (1c
), and generate estimates for the cell densities
and at the points to integrate (1d).The approximation of and at points can be done using a variety of interpolation methods, such as linear or spline interpolation. Since the grid points on which is computed are chosen beforehand to provide a good approximation of , they will generally also form a good set of interpolating points (Franz and Erban, 2013).
The interpolation of the cell densities and from the cell positions to the grid points , necessary to update according to (1d), is slightly more delicate. A basic approach would be to “shift” the delta function from to its closest grid point , so that is a count of the number of type cells in the neighbourhood of . However, this is quite a crude approximation to make, rendering the more accurate approaches in the discretisation of the equations or the interpolation of a waste of effort. This is why the standard approach is to obtain a continuous approximation of the density, or density estimate, and then evaluate it at the grid points . One way to achieve this is to use a particleincell method with piecewise linear polynomials as done by Dallon and Othmer (1997). They use a square lattice and the mass of a delta function at
is distributed among the four nearest grid points proportionally to their distances. Another way is to use a kernel density estimation, which is generally used to estimate the probability density of a random process from a large number of iterations
(for details see Franz and Erban, 2013). The estimate of the density of type cells is found as(7) 
where is the spatial convolution, and denotes a kernel of bandwidth , , taken to be a continuous, symmetric and normalised function. The idea is that, as the bandwidth parameter , the estimate approximates the sum of delta functions in . The Gaussian kernel, , is one of the most commonly used kernels for density estimation. In the context of hybrid modelling, the Gaussian kernel density estimation was used in Franz et al (2013); McLennan et al (2012).
The earliest hybrid models for chemotaxis modelled cells as point particles. Accordingly, the formulation of (1d) with Dirac deltas was considered to be the exact model, and the kernel density estimation its approximation. However, if the actual size of cells is taken into account, it makes sense to assume that the cells consume or degrade chemical all along their shape, and not only in the centre. In this case, the bandwidth parameter can be thought of the lengthscale of the cell (for example, can be related to the hardsphere diameter when cells are modelled as hard bodies), independent of the grid spacing for the chemical (McLennan et al, 2012).
3 Model implementation with Aboria
[basicstyle=,keywordstyle=]‘
We consider an example in a square domain, , with noflux boundary conditions. Cells interact with each other via a soft shortrange repulsive potential , with . Initially there are the same number of particles of each type, , and . Cells of type are distributed according to a twodimensional normal centred at the origin and . Cells of type are uniformly distributed in the whole domain. Initially there is no chemical in the domain.
3.1 Particle set with two types of particles
We define three Aboria variables: ‘type‘ to refer to the cell type (‘type=true‘ for cells of type , and ‘type=false‘ for cells of type ), ‘concentration‘ for the concentration of
at the location of the particle, and ‘drift‘ to store the twodimensional drift vector
(using an Aboria ‘vdouble2‘ type, representing a twodimensional vector). We then define the particle set type, given by ‘Particles_t‘, which contains the Aboria variables and has a spatial dimension of two (specified by the second template argument). For convenience, we also define ‘position‘ as the ‘Particles_t::position‘ subclass (we will use this later on):Finally we create ‘particles‘, an instance of ‘Particles_t‘, containing particles, and initialise the random seed of the set based on a unique sample number ‘sample‘:
In Aboria there are two ways to access and operate on ‘particles‘, either in lowlevel language or a highlevel symbolic language. We show how these two approaches work when initialising the positions of the particles. The lowlevel approach uses the standard C++ random library to generate the gaussian and normal distributions, and loops over the particles to set their positions. We define ‘min‘ and ‘max‘ as vectors representing respectively the lower and upper boundaries for each dimension, ‘min‘ and ‘max‘:
For the highlevel symbolic approach, we define symbolic objects ‘x‘ and ‘typ‘, representing the ‘position‘ and ‘type‘ variables (and others that we will use later on). We also define two Aboria random number generators; ‘normal‘ for normally distributed numbers, and ‘uni‘ for uniformly distributed numbers. Finally, we also create a label object ‘k‘ associated to the particle set. This label performs a similar function to the subscript for the variable in Eq. 2. All operations involving ‘k‘ are implicitly performed over the entire particle set.
Then the highlevel initialisation of the positions of ‘particles‘ (assuming that the ‘type‘ variable has already been set) is:
The advantage of the highlevel approach is that we can directly write expressions which are meant to be evaluated over the whole particle set. However, as it implements a custom eDSL, it is by definition limited to operations that can be expressed by this language. For example, the interaction between the chemical grid and the individual particles cannot be expressed with this DSL. Therefore, for the remainder of this paper we will proceed by implementing the model using the lower level interface.
Finally, we initialise the spatial search data structure, providing it with lower (‘min‘) and upper (‘max‘) bounds for the domain, and setting nonperiodic boundary conditions (‘periodic‘ in this case is ‘false‘):
This subdivides the computational domain into square cells of equal size. The default cell side length is such that, if particles were uniformly distributed in the domain, then there would be on average ten particles per cell. However, it is also possible to specify a desired side length as a forth parameter (to make it equal, for example, to a cutoff distance for the calculation of interaction forces).
3.2 Equation of motion of cells
Here we explain how to implement the SDEs (1a) and (1b). For simplicity, in this section we are going to assume that the chemical drift in (1b) is fixed and already stored in the ‘drift‘ variable within each particle. To implement the interaction force between two particles and , we use the ‘euclidean_search‘ function in Aboria, which returns an iterator ‘j‘ that iterates through all the particles within a certain radius (‘cutoff‘) of a given point. For each potential pair of particles, the iterator ‘j‘ can provide the shortest vector between and (normally , but not always for periodic simulations), and we use this vector to evaluate the interaction force between and :
Note that we increment the variable ‘next_position‘, rather than ‘position‘, as this interaction loop uses the particle positions and therefore cannot update them until the loop is complete.
Now we can implement both the drift and the Brownian diffusion of a particle using the stored ‘drift‘ variable, as well as a random number generator that is stored in the variable ‘generator‘. We will use the standard C++ normal distribution (the ‘normal‘ variable) to generate two normally distributed random numbers used for the diffusion:
The noflux boundary conditions are enforced by reflections if the particles end up outside the domain:
3.3 Hardsphere interactions
In the previous subsection, we showed how to implement soft interactions between particles. Below we show an example of the update of ‘x‘ if, instead of using an interaction potential, we want to model cells as hard spheres of diameter ‘epsilon‘. Using the first of the two hardsphere collision algorithms discussed, the new particle positions are obtained using:
3.4 Reactions between cells
The reactions (1c) to change type between cells are implemented as follows:
If the ‘type‘ variable is ‘true‘ ( cell), then the cell changes type if , where . If instead ‘type‘ is ‘false‘ ( cell), then the reaction takes place if . In our implementation, the chemical concentration at the position of the th particle is saved in the ‘conc‘ variable (we explain in the next section how this is updated).
3.5 Chemical concentration field
The chemical is modelled by its continuum concentration rather than individual particles. For this reason, instead of using an Aboria ‘Particles_t‘ set to describe it, we use a standard C++ linear algebra library. We take the computational domain for the chemical to be equal to that of the cells, that is, . We discretise the domain in grid points, , in each direction, and integrate the PDE (1d) for the chemical using finite differences (secondorder centred differences in space, and forward Euler in time). ‘Vector‘ is a vector type of length used to create instances ‘c‘ for the discretised chemical concentration , ‘rhoalpha‘ for the type cell density estimate , and ‘rhobeta‘ for the type cell density estimate . ‘Vector2‘ is a matrix type of size to store the gradient of the chemical concentration . Finally, ‘SparseMatrix‘ is a sparse matrix type of size that is used to store the various finite differences discretisation matrices: ‘D2‘ stores , where is the discrete Laplacian and is the timestep, and ‘D1x‘ and ‘D1y‘ that perform the first derivatives with respect to the horizontal and vertical coordinates, respectively:
Then we define the Gaussian kernel to estimate the cell densities and , rescaled with (the cell diameter) since we assume cells produce/consume chemical using receptors that are distributed over their bodies (McLennan et al, 2012):
We approximate the chemical concentration and the cell density estimates at the grid points. The position of the grid point closest to the left bottom corner of the computational domain is and, as mentioned above, is stored as ‘min‘. The code below computes the discretised cell density estimates, finding the cells of a given type that are within a ‘cutoff‘ distance of the grid point where we want to compute ‘rhoalpha‘ and ‘rhobeta‘ at. Then we evaluate the kernel at and add it to the corresponding vector and grid position:
We implement (1d) to update the chemical concentration vector ‘c‘ using:
Then we calculate the concentration gradient ‘grad_c‘ like so:
3.6 Spatial matching from regular grid for the chemical to particle positions
Finally, we need to convert the continuum variables and (which in the numerical simulation are approximated at regular grid points) to the positions of type cells, so that we can evaluate the drift in (1b) and the reaction rate in (1c). To do that, we use linear interpolation:
4 Simulation results: output and analysis of data
While the C++ language is ideal for implementing fast simulations with low memory overhead, the Python language is generally preferred for plotting, and pre and postprocessing of simulation data. Thankfully, there are many tools that enable wrapping of C++ code in Python, and we will make use of one of these, Boost Python, to enable us to call our simulation code from Python.
The main difficulty in wrapping C++ code in Python is transferring data between the two languages. In our case, we only need to transfer data to Python for postprocessing and plotting. We can use the Boost NumPy extension to wrap a Numpy array around a Aboria variable given in the template argument ‘V‘.
Note that no copying of data occurs in this function, the new Numpy array that is returned from the function simply wraps the data so that it can be easily accessed in Python. This function assumes that the Aboria variable is a vector type (e.g. ‘position‘), but we can easily write another function to wrap a scalar variable (e.g. ‘type‘). Note also that we can write a very similar function to transfer the grid data to Python as well (see the full code in the Supplementary Material for all three of these functions).
Now that we can transfer data, we need a C++ object with which we can interact in Python. Thus we will create a ‘Simulation‘ class to store our data, with functions like ‘integrate()‘ and ‘get_positions()‘ that will allow us to either integrate the simulation forward in time, or obtain internal variables for plotting:
Now we can use Boost Python to wrap our ‘Simulation‘ class and enable it to be used from Python:
After this we need to compile the code we gave generated thus far. For the example code included with this paper we have used the CMake build system (see the ‘CMakeLists.txt‘ file included in the Supplementary Material for details of how this is done). After compilation we are left with a final library file named ‘chemo‘ (with an extension that depends on the particular operating systems you are using) that we can use in Python like so:
The above listing simple integrates the simulation until , and then creates a scatter plot of the cells coloured by their type. This simulation consists of a single sample, and creating new simulation objects with a different initial seed (in the code above the sample seed is set to 1) will result in a different random realisation.
In order to gain a complete picture of the dynamics we need to run multiple samples, and average the results. However, now that our simulation code is running in Python, we can use its highlevel features and libraries to do this relatively easily. For example, we can use Python’s standard multiprocessing library to run many simulations in parallel, average the particle positions by calculating histograms using Numpy’s histogram2d function, cache the results of the simulations to disk using Python pickle, and finally plot the results using matplotlib. We will not explain these facilities in detail here, but instead refer readers to the external documentation links provided, and to the Supplementary Material for a full code listing in ‘paper_plots.py‘ showing how this might be done.
The generated figures from ‘paper_plots.py‘ are plotted in Figure 1, showing the averaged histograms for the and particle types, as well as the averaged distribution of the chemical concentration . The upper row shows the diffusion of particles from an initial Gaussian profile becoming increasingly spread out over the domain. The middle row shows the diffusion of particles from an initial uniform gradient. As the simulation proceeds the particles produce the chemical in the centre of the domain, and soon the chemotactic gradient term results in particles becoming clustered in the middle of the domain. The influence of the noflux boundary conditions is seen as a buildup of particles near all four boundaries.
Figure 2 shows the number of and particles, as well as the total number of all particles. This shows the dominant conversion of particles to , as well as the conservation of over time.
Once we have the simulation framework established, we can begin conducting numerical experiments, altering the domain, parameters, initial distributions or behaviours of particles in order to explore different chemotactic behaviours, or taking advantage of the transparency of numerical simulation by measuring different quantities of interest. For example, in order to explore the average mean squared displacement for the and particles we might wish to run another experiment using periodic boundary conditions, and track the total displacement of an individual particle from its initial position. We can do this by defining a new variable ‘starting‘:
and then initialising it to the starting position of each particle at the beginning of the simulation.
Instead of the noflux boundary conditions that we implemented previously (in Section 3.2), we will use periodic boundary conditions. Note that setting the periodic argument of ‘init_neighbour_search‘ to ‘true‘ (see Section 3.1) will cause particles that cross the periodic boundary to be automatically moved to the opposite boundary. We will counteract this by updating ‘starting‘ whenever a particle will cross the periodic boundary, so that the particle displacement thus far is preserved.
Once we have our new variable in place we can simply calculate the mean squared displacement (MSD) of each particle by comparing its current ‘position‘ to the ‘starting‘ variable. We can use our Python wrapper to obtain both of these variables and calculate the MSD as a postprocessing step.
Figure 3 shows the MSD for three different scenarios: (1) , (2) , and (3) . In all three the chemical gradient was set to be a constant by setting the chemical diffusion, production and consumption to zero (), and the initial profile equal to . The chemical concentration does not satisfy the periodic boundary conditions in this case, however, its gradient is does and this is the only factor influencing the simulation via the drift term on the particles. The diffusion constant for each species was , and reactions between particles types are turned off. For simulation (1) the particles are not affected by the drift term and so the MSD is a straight line with a low gradient. For simulation (2) the particles are affected by the constant drift, which adds a dominant term to their MSD. For simulation (3) the reactions are turned on again and there are both and cells, so the net effect on the MSD is a combination of these two behaviours.
5 Conclusions
Particlebased models for biological processes have become of widespread use in mathematical biology. These come in many forms, with particle motions described by discrete or continuum random walks, complex interactions between particles (including reactions, hardcore interactions, or interaction potentials), and interactions with the environment (such as chemotaxis or transport through crowded or heterogeneous domains). In some cases, one also requires to couple a particlebased model, describing for example individuals cells, with a PDE model to represent a continuum field (such as a chemical concentration). In addition, despite their apparent simplicity, particlebased models can be challenging to implement and simulate, as they tend to scale badly with the number of particles in the system (which can be large in many applications) and, due to stochasticity, often many realisations of the same simulation are required.
This diversity in particlebased models combined with the computational challenges in simulations, makes the implementation of particlebased models far from straightforward. In this paper we presented Aboria, a C++ library, designed to provide the flexibility required to implement particlebased models commonly used in mathematical biology, in a high performance and easy to use fashion.
We have demonstrated the usage of Aboria implementing a model for cell diffusion and chemotaxis with shortrange interactions. The model has many of the features described above, namely, cells move according to biased Brownian motion, they interact with each other and with a chemical (that is modelled as a continuum), and there are reactions that change the number of particles. We have shown how Aboria can be used in combination with Python to produce outputs such as the cell densities, numbers, and the mean square displacement.
References
 Abraham et al (2015) Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindahl E (2015) Gromacs: High performance molecular simulations through multilevel parallelism from laptops to supercomputers. SoftwareX 1:1925
 Andrews and Bray (2004) Andrews SS, Bray D (2004) Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys Biol 1(3):137151
 Betzig et al (2006) Betzig E, Patterson GH, Sougrat R, Lindwasser OW, Olenych S, Bonifacino JS, Davidson MW, LippincottSchwartz J, Hess HF (2006) Imaging intracellular fluorescent proteins at nanometer resolution. Science 313(5793):16421645
 Bruna (2012) Bruna M (2012) Excludedvolume Effects in Stochastic Models of Diffusion. PhD thesis, University of Oxford
 Bruna and Chapman (2012) Bruna M, Chapman SJ (2012) Excludedvolume effects in the diffusion of hard spheres. Phys Rev E 85(1):011103
 Bruna et al (2017) Bruna M, Chapman SJ, Robinson M (2017) Diffusion of Particles with ShortRange Interactions. SIAM J Appl Math 77(6):22942316
 Dallon and Othmer (1997) Dallon JC, Othmer HG (1997) A discrete cell model with adaptive signalling for aggregation of Dictyostelium discoideum. Phil Trans R Soc B: Biol Sci 352(1351):391417
 D’Orsogna et al (2006) D’Orsogna MR, Chuang Y, Bertozzi AL, Chayes LS (2006) SelfPropelled Particles with SoftCore Interactions: Patterns, Stability, and Collapse. Phys Rev Lett 96(10):104302
 Eastman et al (2012) Eastman P, Friedrichs MS, Chodera JD, Radmer RJ, Bruns CM, Ku JP, Beauchamp KA, Lane TJ, Wang LP, Shukla D, et al (2012) Openmm 4: a reusable, extensible, hardware independent library for high performance molecular simulation. Journal of Chemical Theory and Computation 9(1):461469
 Erban et al (2007) Erban R, Chapman J, Maini PK (2007) A practical guide to stochastic simulations of reactiondiffusion processes. arXiv 0704.1908
 Fong and Darve (2009) Fong W, Darve E (2009) The blackbox fast multipole method. Journal of Computational Physics 228(23):87128725
 Franz and Erban (2013) Franz B, Erban R (2013) Hybrid Modelling of Individual Movement and Collective Behaviour. In: Dispersal, Individual Movement and Spatial Ecology, Springer, Berlin, Heidelberg, Berlin, Heidelberg, pp 129157
 Franz et al (2013) Franz B, Xue C, Painter KJ, Erban R (2013) Travelling Waves in Hybrid Chemotaxis Models. Bull Math Biol 76(2):377400
 Gillespie (1977) Gillespie DT (1977) Exact stochastic simulation of coupled chemical reactions. J Phys Chem 81(25):23402361
 Greengard and Rokhlin (1987) Greengard L, Rokhlin V (1987) A fast algorithm for particle simulations. Journal of computational physics 73(2):325348
 Guo et al (2008) Guo Z, Sloot PMA, Tay JC (2008) A hybrid agentbased approach for modeling microbiological systems. J Theor Biol 255(2):163175
 Hutzenthaler et al (2012) Hutzenthaler M, Jentzen A, Kloeden PE (2012) Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients. The Annals of Applied Probability 22(4):16111641
 Jeon et al (2010) Jeon J, Quaranta V, Cummings PT (2010) An OffLattice Hybrid DiscreteContinuum Model of Tumor Growth and Invasion. Biophys J 98(1):3747
 McLennan et al (2012) McLennan R, Dyson L, Prather KW, Morrison JA, Baker RE, Maini PK, Kulesa PM (2012) Multiscale mechanisms of cell migration during development: theory and experiment. Development 139(16):29352944
 Middleton et al (2014) Middleton AM, Fleck C, Grima R (2014) A continuum approximation to an offlattice individualcell based model of cell migration and adhesion. J Theor Biol 359:220232
 Newman and Grima (2004) Newman TJ, Grima R (2004) Manybody theory of chemotactic cellcell interactions. Phys Rev E 70(5):051916
 Robinson and Bruna (2017) Robinson M, Bruna M (2017) Particlebased and meshless methods with Aboria. SoftwareX 6 IS :172178
 Scala et al (2007) Scala A, Voigtmann T, De Michele C (2007) Eventdriven Brownian dynamics for hard spheres. J Chem Phys 126(13):134109
 Stroustrup (2013) Stroustrup B (2013) The C++ programming language. Pearson Education
Comments
There are no comments yet.