The auxiliary region method: a hybrid method for coupling PDE- and Brownian-based dynamics for reaction–diffusion systems

Reaction–diffusion systems are used to represent many biological and physical phenomena. They model the random motion of particles (diffusion) and interactions between them (reactions). Such systems can be modelled at multiple scales with varying degrees of accuracy and computational efficiency. When representing genuinely multiscale phenomena, fine-scale models can be prohibitively expensive, whereas coarser models, although cheaper, often lack sufficient detail to accurately represent the phenomenon at hand. Spatial hybrid methods couple two or more of these representations in order to improve efficiency without compromising accuracy. In this paper, we present a novel spatial hybrid method, which we call the auxiliary region method (ARM), which couples PDE- and Brownian-based representations of reaction–diffusion systems. Numerical PDE solutions on one side of an interface are coupled to Brownian-based dynamics on the other side using compartment-based ‘auxiliary regions’. We demonstrate that the hybrid method is able to simulate reaction–diffusion dynamics for a number of different test problems with high accuracy. Furthermore, we undertake error analysis on the ARM which demonstrates that it is robust to changes in the free parameters in the model, where previous coupling algorithms are not. In particular, we envisage that the method will be applicable for a wide range of spatial multi-scales problems including filopodial dynamics, intracellular signalling, embryogenesis and travelling wave phenomena.


Introduction
Reaction-diffusion models are important mathematical tools that are used to represent and understand complex biological and physical behaviours. They model the random movement of the particles (diffusion) and the interactions between particles (reactions), giving them a wide array of applications across multiple spatial scales. These applications range from the large-scale representation of striped vegetation patterns in semi-arid landscapes [1] and the spread of epidemics [2] to smaller-scale studies of pattern formation during embryogenesis [3,4] and, at even smaller scales, to the study of actin dynamics inside a cell's filopodia [5] and intracellular dynamics [6][7][8].
Reaction-diffusion models can be specified at different levels of detail depending on the temporal, spatial and concentration scales involved in the application (table 1). At the finest scale that we will consider are microscopic dynamics. These models and methods (which include Brownian motion for purely diffusive systems and Smoluchowski dynamics [9,10] or Green's function reaction dynamics (GFRD) [11] for reaction-diffusion systems) are among the more detailed representations of such systems, but consequently are relatively computationally expensive. 1 They require not only the knowledge of the location of all particles at all times, but in the case of second-and higher-order reactions, the pairwise distances between particles, which requires large memory, and are expensive to calculate for many time-steps. In the case of diffusion-limited reactions, time-steps must be taken to be extremely small to ensure that reactive particles do not jump past each other and that the attendant reaction events are not missed. All update steps also require the production of a normally distributed random number for each coordinate of each particle which can be computationally expensive depending on the reaction system that is being modelled. However, some of these expensive steps can be accelerated by considering event-driven algorithms or employing approximate algorithms with longer time-steps. GFRD is an event-driven algorithm that differs from the standard method for simulating Brownian motion. It uses a maximum time-step so that only single particles, or pairs of particles, need to be considered. It then uses the exact solution to the Smoluchowski equation in order to combine the movement of, and interactions between, particles. If particles are far apart, the event-based time-steps are large. Smoldyn uses relatively long time-steps, and accounts for the error that this causes (due to possible reactant pairs passing by one another without the possibility of reacting) by making the effective particle sizes larger. Microscale modelling is particularly useful when fine scale detail is required, for example, when considering the binding of particles to receptors [12][13][14]. An even finer scale representation, in which atomistic dynamics can be represented, is available, if required. Typically, modelling at this scale is known as molecular dynamics, and we direct the reader to Holley [15] and Dürr et al. [16] for more information about modelling at this scale.
At a coarser scale, we have compartment-based or mesoscopic models. Like the fine-scale, microscopic models, these also account for stochastic variation; however, particles are now considered to belong to compartments rather than having their exact locations tracked. Particles can either react with one another within a compartment, or can jump between adjacent compartments with given rates, simulating diffusion. Compartment-based models can be simulated using either exact but computationally expensive [17][18][19] or inexact but computationally cheaper [20] stochastic simulation algorithms (SSAs). The exact methods (so-called because they produce sample paths consistent with the associated chemical master equation (CME)) effectively assign exponential waiting times to every possible event (diffusive jump or reaction) and then choose the event with the shortest waiting time to enact. In general, they are faster than the microscopic methods, as pairwise reaction distances do not need to be calculated for bimolecular reactions and individual particle identities are not tracked, but are less accurate, because they only record a particle's location up to the accuracy of the compartment size, and generally particles are only allowed to react with others in the same compartment [21].
Finally, at the coarsest scale lie continuum or macroscopic models. The most commonly employed macroscopic models for reaction-diffusion systems comprise partial differential equations (PDEs). 2 These methods are generally only valid for high particle numbers. The stochastic variations, which are considered small enough to be neglected at high copy numbers, play a pivotal role in the dynamics at low copy numbers, leading the PDE solutions to diverge from the true underlying dynamics. There is a wealth of well-established numerical methods that can quickly simulate an approximate solution to a  PDE. These include finite-difference methods, finite-volume methods and finite-element methods (e.g. [22][23][24][25]). Often though, important biological and physical phenomena are genuinely multiscale [26][27][28][29]. In spatial reaction-diffusion systems, concentration may vary over orders of magnitude. In regions of low concentration, it is often important to employ detailed individual-based models in order to correctly represent the dynamics. If these models were to be employed indiscriminately throughout the domain, however, the regions of high concentration, in which there are many individual particles to be evolved, might render the system computationally intractable. In these regions it might be acceptable to employ a coarser and less computationally expensive model. A canonical example of this phenomenon is the stochastic Fisher wave [30,31]. The wave speed is determined by the stochastic activity at the pulled front, so it is important to employ an accurate individual-based representation of the dynamics in this region. Conversely, behind the wavefront, the detailed dynamics are of little importance. It is possible, therefore, to employ a coarser, cheaper representation of the dynamics in this region.
Spatially coupled hybrid methods have been developed for precisely this purpose: to simulate spatially inhomogeneous domains both accurately and efficiently. In general, such methods are designed to accelerate expensive computations whilst maintaining reasonable levels of accuracy. The majority of spatially coupled hybrid methods divide the computational domain into distinct regions using interfaces. The dynamics of adjacent regions are represented using different methods. Regions in which detailed representations of the dynamics are required for accuracy are simulated using a fine-scale method, whereas regions in which less detail is required are modelled with a coarser, less computationally expensive method. There can be two reasons for this. The first is in order to resolve a particular region of the spatial domain in more detail, such as when looking at the behaviour of ions around gated channels [12], or when building a model for the energy in a liquid crystal [32]. Both of these examples have a prohibitively slow but accurate model that is required in certain regions of space, but which is too computationally expensive if used everywhere. The second reason is to simply segregate a region of the domain in which there are very few particles. In these regions, a coarse method (for example a continuum model) may be too inaccurate.
There exist hybrid methods that couple each of the different scales described above to one another (and indeed many more; see Smith & Yates [33] for a comprehensive review of such methods). Macroscopic-to-mesoscopic methods have been proposed which employ averaged fluxes in order to calculate appropriate boundary conditions for each regime at the interface(s) [34][35][36], as well as using an extra compartment within the macroscopic region [37]. Mesoscopic-to-microscopic methods, which also employ extra compartments, this time in the microscopic regime, have been developed [38]. Moreover a class of methods using adapted rates of diffusion from the mesoscopic to the microscopic domains have been proposed and successfully applied to represent biological processes [5,12,29,39,40]. There are fewer macroscopic-to-microscopic hybrid methods in the literature. Macro-to-micro methods that allow mass to flow over the interface in both directions in order to initialize particles [41] or that average solutions on either side of the interface to find a flux [42] can be found in the literature. For a more detailed review of spatially extended hybrid methods, see [33].
Two of the above-mentioned hybrid methods are of particular relevance for the purposes of this paper. The pseudo-compartment method, presented by Yates & Flegg [37], is a macroscopic-tomesoscopic (specifically PDE-to-compartment) method in which the coupling is achieved using an extra compartment, known as the 'pseudo-compartment', adjacent to the interface within the macroscopic  Figure 1. Schematics for (a) the pseudo-compartment method and the (b) the ghost cell method. The green line represents the PDE solution, the blue boxes the number of particles in the mesoscopic region of the respective hybrid methods, and the yellow dots denote the Brownian particles. These particles are shown with a volume, but in the simulations do not have a mass or volume. In the scenario illustrated, the particles reside on the one-dimensional line, but have been illustrated in the plane in order to show the directions and magnitudes of their next movement clearly (black arrows). The green boxes in (a) denote the number of particles in the pseudocompartment, and similarly, the yellow boxes in (b) are the number of particles in the ghost cell, with each box representing a single particle. In each case, the red line denotes the point interface between the two regimes.
domain. In this compartment, mass is represented using both the PDE solution and the compartmentbased method (with particle numbers found by direct integration of the PDE over this region). Particles are then allowed to cross the interface in both directions using the compartment-based method. We give a schematic of this method in figure 1a.
The ghost cell method proposed by Flegg et al. [38] is a mesoscopic-to-microscopic method which uses an extra compartment in the microscopic domain. The number of particles in this 'ghost cell' is simply the number of Brownian particles which reside in this region. Again, particles are allowed to jump across the interface using the compartment-based mesoscopic method. A schematic of the method is given in figure 1b.
In this paper, we employ the two methods described above (see figure 1) in order to couple a macroscopic PDE description for reaction-diffusion systems to a corresponding microscopic Brownian dynamics representation through the use of 'auxiliary regions'. These regions are compartments, which lie either side of the interface, and allow mass to pass between the two regimes via a mesoscopic jump process (see figure 4 for a schematic). Within the auxiliary regions, mass is simultaneously represented using both the description for the region in which they reside (i.e. PDE or Brownian) and the mesoscopic description. Changes (i.e. reactions or diffusion events) implemented under one modelling paradigm (e.g. the compartment-based representation of the auxiliary region) are simultaneously implemented in the other (e.g. the PDE or Brownian representations in these regions). The interface, which divides the two modelling paradigms, can either be static, in which case it remains in its initial position, or adaptive, in which case it moves with the density profile in order to ensure that regions of space with few particles are simulated using the finest scale. Through a series of test cases, we demonstrate our algorithm to be more accurate and more robust to model parameters than previous PDE-to-Brownian coupling algorithms.
The paper is organized as follows. In §2, a previous attempt at hybridizing a Brownian dynamics model to its corresponding mean-field PDE description is evaluated in more detail [41]. A description of our novel auxiliary region method (ARM) is presented in §3 alongside the relevant justifications and pseudocode. Numerical results, verifying the accuracy of our hybrid method, are presented in §4. Numerical error analysis is conducted in §5, where we also discuss restrictions on the model parameters for the effective functioning of the coupling algorithm. We conclude with a discussion of the effectiveness of our new hybrid method and suggest avenues for further exploration in §6.

An existing PDE-to-Brownian coupling
In this section, we summarize the pioneering work of Franz et al. [41], who were among the first to couple PDE and Brownian dynamics representations of reaction-diffusion. By replicating their results,    [41]. Descriptions are as in figure 1. The PDE mass labelled α (orange) is the density of PDE mass that has flowed over the interface in the Brownian update step. The peak in the PDE curve near the interface represents the addition of a δ-function corresponding to a Brownian particle that crosses the interface.
we demonstrate that their 'PDE-assisted Brownian dynamics' algorithm is not robust to simulation parameter choice, even for simple diffusive processes. This motivates the need for a more robust coupling method, which we provide in the form of the ARM in §3.

PDE-assisted Brownian dynamics
Hybrid methods that couple the PDE description of a reaction-diffusion system to its corresponding Brownian dynamics representation have been relatively poorly investigated in comparison to PDE-tocompartment-based and compartment-based-to-Brownian couplings. In part, this is a result of the fact that such hybrid algorithms neglect mesoscale representations of particle dynamics, meaning that they must bridge a greater scale separation than either of the other two hybrid paradigms. Mainly though, the absence of many examples of PDE-to-Brownian hybrid methods is due to the inherent difficulty when converting PDE mass to individual particles (and vice versa) when coupling Brownian dynamics models to continuum PDE representations. Below, we describe two algorithms proposed by Franz et al. [41], but focus on the first, a method with an interfacial coupling. We choose to focus on this coupling because our ARM coupling method, described in §3, also uses an interface. Franz et al. [41] present two related algorithms. In the first, the non-overlapping PDE and Brownian domains are separated by an interface (figure 2). Both PDE and Brownian representations are updated using a time-driven algorithm, with the PDE time-step much smaller than the Brownian time-step. The discretized PDE is evolved (until the time reaches the next Brownian time-step) using a centred finitedifference scheme with implicit Euler time-stepping, and PDE mass is allowed to cross the interface between the two regimes. Provided that the Brownian time-step is sufficiently small, the amount of mass that crosses the interface between Brownian time-steps gives the probability that a new particle is placed within the Brownian domain. A uniformly distributed random number is used to determine whether a particle is initialized in the Brownian regime or not. If it is, this particle's position is randomly initialized according to the normalized density profile of the PDE mass that crossed the interface in the previous Brownian time-step. If a Brownian particle crosses into the PDE domain, a particle's worth of mass is added to the PDE solution at its new location as a δ-function and the individual particle is removed. We have illustrated this method schematically in figure 2.
Franz et al. [41] found the variance in particle numbers in the Brownian region of the hybrid domain to be altered in comparison to the variance that would be expected in a fully Brownian simulation. In order to counteract this problem, they introduced a second algorithm, in which an overlap region replaces the interface. Within the overlap region, mass can be simulated as either Brownian particles or as part of the PDE. The coupling works in the same way as in the interfacing algorithm; however, the Brownian particles are subsumed into the PDE only once they have crossed the boundary of the overlap region closest to the fully PDE domain. Similarly, PDE mass can only be converted to Brownian particles once it has flowed over the overlap boundary adjacent to the fully Brownian domain.
The Brownian time-step in the algorithm is required to be small, in order that the total probability of initializing a particle in the Brownian regime is less than one. However, the algorithm runs into  difficulties if the time-step is chosen to be too small. Specifically, the amount of mass that flows over the interface between updates of the Brownian dynamics is too small in comparison to that which would be predicted theoretically using the exact diffusion kernel. This gives rise to inaccuracies in the algorithm, particularly if long simulation times are required. This sensitivity to the choice of Brownian time-step restricts the physical scenarios to which the algorithm can be applied.
In figure 3, we present three snapshots of the evolution of the first version of the algorithm (interface rather than overlap region) which illustrate this problem. By time t = 2, in figure 3c, there is a clear disparity between the hybrid method and the mean field solution (black dotted line). Disparities of this nature are not acceptable when modelling real reaction-diffusion systems, irrespective of the computational savings the algorithm is able to produce.

The auxiliary region method
In this section, we present our novel ARM for coupling PDE and Brownian-based representations of reaction-diffusion. For simplicity, we will present a version of the method with a single interface separating two regimes. However, the method can be easily generalized to multiple interfaces which separate alternating PDE and Brownian regions. Sequentially, we describe the composition of the domain and the models we employ in each region; the nature of the auxiliary regions; the implementation of movement of mass across the boundary; the implementation of reactions; and finally the specific details required for the simulation of the algorithm, including pseudocode for its implementation. All code, which has been written in Matlab, can be found in the electronic supplementary material.   figure 1. The interface is the red line in the centre and the two auxiliary regions are shown with blocks to indicate the number of particles residing within them. In the PDE and Brownian auxiliary regions, each block represents a particle in the compartment-based representation and the number of blocks is determined by integrating the PDE over the auxiliary region Ω PA , and counting the number of Brownian particles in Ω BA , respectively.

The domain composition
Recall that, for our coupling method, space is partitioned into two regions within which we use different modelling paradigms (PDE and Brownian dynamics) to simulate the underlying reaction-diffusion system. Separating the two regions is a point interface, over which particles can jump according to a compartment-based method.
Consider a one-dimensional domain 3 We split Ω into two regions, Ω P = (L 1 , 0) and Ω B = (0, L 2 ) (separated by an interface I at position 0), within which the evolution of the system will be represented using a PDE description and Brownian dynamics, respectively.

The auxiliary regions
Particles can move between the two domains (Ω P and Ω B ) via the auxiliary regions Ω PA and Ω BA ; subsets of Ω P and Ω B respectively, each of width h a > 0. Within these regions, mass/particles are simultaneously represented according to the default methodology for their domain (either PDE in Ω P or Brownian dynamics in Ω B ), but also as well-mixed particles in their respective auxiliary regions Ω PA and Ω BA . These auxiliary regions act as a bridge between the fine-and coarse-scale descriptions. A schematic of the domain's composition is given in figure 4.
We justify the use of the Brownian auxiliary region by following the methodology set out in Flegg et al. [38]. The entire Brownian domain can be simulated using a mesoscopic compartment-based regime, and equivalently using a microscopic simulation. In the absence of reactions, if the particles in the microscopic simulation are 'binned' into the same compartments as the mesoscopic simulation, the expected numbers in each compartment for each simulation would be the same. At this scale, the two methods are equivalent ways of simulating the same diffusive process [38].
To justify the use of the PDE auxiliary region, we appeal to the arguments of Yates & Flegg [37]. We note that the PDE density can be thought of as the probability of finding a particle at a particular position and time, scaled by the number of particles within the PDE domain. Provided that the auxiliary region is sufficiently narrow, the PDE density within the auxiliary region can be thought of as being approximately uniformly distributed across the region with the appropriate number of particles. This is precisely the interpretation of the contents of a compartment within the mesoscopic, compartment-based framework.

The PDE regime, Ω P
Within Ω P , we represent the mass of particles using and t)) T , denotes the density of species k = 1, . . . , K at position x and time t, D is a diagonal matrix containing the Fickian diffusion constants for each species, and f is a function that encapsulates the effect of any reactions on each species. We also use the notation ∂Ω P to represent the boundary of Ω P , and c 0 (x) is the initial condition. For all the simulations presented in this paper, we employ the finite-difference θ -method (a general family of finite-difference methods). 4 Although the Crank-Nicolson method (θ = 0.5) is second-order accurate and unconditionally stable, we use θ = 0.51 because the Crank-Nicolson method can give rise to spurious oscillations when implemented on step-function initial conditions of the sort we will consider [22].

The Brownian regime, Ω B
Within Ω B , all particles are tracked and their positions updated according to the following (computational) stochastic differential equation which simulates Brownian motion: where y k i (t) denotes the location of particle i of species k within Ω B , t is the time-step for both the PDE and Brownian dynamics simulators 5 and N k HB (t) is the number of particles of species k in Ω B at time t. Once again, we set reflective boundary conditions at both ends of Ω B to ensure that no particles can leave this domain via a Brownian diffusion event. The zero-flux boundary conditions at the interface for both PDE and Brownian regimes ensure that mass can only cross the interface according to the compartmentbased method.

Movement across the interface
As both domains, Ω P and Ω B , have zero-flux boundaries at the interface, particles can only cross over the interface via the auxiliary regions. In effect, these regions comprise a two-compartment reaction-diffusion master equation (RDME) model. Each particle in each auxiliary region jumps to its neighbouring region on the other side of the interface with a rate d k (for species k), which is related to the macroscopic diffusion coefficient (for species k), D k , via Here, h a is the width of each auxiliary region, which is assumed to be the same for both the Brownian and PDE auxiliary regions. In order to implement jumps (or reactions, where necessary) according to the RDME, we require particle numbers. Borrowing terminology from Yates & Flegg [37], the number of 'pseudo-particles' of species k within the PDE auxiliary region, Ω PA , at time t, denoted N k PA (t), is calculated as The number of particles of species k in the Brownian auxiliary region, Ω BA , is given by These particle numbers allow us to define propensity functions corresponding to diffusive jumps between, or reactions within, the auxiliary regions. For diffusive jumps between the two auxiliary regions, the propensity functions for species k within the PDE and Brownian auxiliary regions are (respectively) α k We note here that if N k PA (t) < 1, we set α k P (t) = 0 to prevent the possibility of negative density. While it may be a problem if this scenario occurs persistently, practically speaking, we should choose the position of the interface such that density is always large enough that this does not happen. An adaptive interface will allow us to satisfy this criterion (see §4.4.3), and hence this problem would not occur when using such an interface.
When a particle jumps from Ω BA to Ω PA , a particle within the Brownian auxiliary region is chosen uniformly at random to be removed, and a particle's worth of mass is added to the PDE solution uniformly across Ω PA for the species, k, which has changed: Similarly, if a jump is enacted in the opposite direction, from Ω PA to Ω BA , we first remove a particle's worth of mass uniformly from Ω PA for the appropriate species k: and a new particle is initialized within the Brownian auxiliary region, Ω BA , with position chosen uniformly at random.

Reaction implementation
Throughout Ω P , all reactions are implemented using the reaction operator f (c). The method we employ to implement reactions within Ω B depends on the location of the reactant particles. Let R denote the set of reaction pathways (with |R| = R). Define the subset of reactions R * (t) at time t as follows: R * (t) = {all reactions for which at least one set of reactant particles lies exclusively within Ω BA }.
Reactions between molecules for which at least one of the reactive molecules lies within Ω B \Ω BA are implemented using an appropriate microscopic approach, such as the λ-ρ method [43,44]. However, if at least one set of participating particles lie in Ω BA (i.e. r ∈ R * ), care needs to be taken over the interaction of such particles and the mass on the other side of the interface in Ω P . As explained below, we will implement the reactions r ∈ R * for these reactant particles using the compartment-based method.
For illustrative purposes, consider a reversible second-order reaction involving species A, B and C: Under the λ-ρ method [43] and its later modification [44], for the forward reaction, a particle of species A and a particle of species B are required to be within a distance ρ of one another in order to react. They then react with a rate λ, where λ is a function of both the reaction radius ρ and the reaction rate κ 1 . Imagine that an A particle (without loss of generality) in Ω B is close enough to the interface that the reaction radius ρ is larger than the distance between itself and the interface. For consistency with the Brownian representation, the A particle should be allowed to react with a B particle in the PDE region.
The implementation of such reactions would be extremely difficult. Instead, by ensuring bimolecular reactions within the auxiliary region are implemented according to the mesoscopic compartment-based method, we avoid such issues (provided that the width of the auxiliary region is chosen to be larger than the interaction radius ρ).
According to the backwards reaction, two particles are created after the reaction has occurred. These particles are placed a certain distance away from each other (called the dissociation radius) in order to achieve a specified probability of geminate recombination (a recombination of any pair of A and B particles that were initialized from the same C particle). If this radius intersects with the PDE regime, then there is the potential for individual particles to be initialized within Ω P . By again employing the mesoscopic representations for reactions, we resolve this issue. All product particles are assumed to be placed uniformly throughout the Brownian auxiliary region. Particles that are products of the backward dissociation reaction in Ω B \Ω BA are extremely unlikely to be placed in Ω P (again, providing that the auxiliary region is larger than the dissociation radius).
For these reasons, all of the reactions r ∈ R * (for which at least one set of participating particles lie in Ω BA ) are implemented using the compartment-based method, in which reactions are incorporated as events in the associated Markov chain, according to the RDME. We can write the following propensity functions for reactions within Ω BA : α r (t) = g r (N BA (t))κ r h 1−ν a , (3.11) for any reaction channel r ∈ R * (t) of order ν and corresponding reaction rate κ r , where N BA (t) = (N 1 BA (t), . . . , N K BA (t)) T and g r is the appropriate number of possible combinations of the reactants for reaction r from the particles that lie within Ω BA . Recall, however, that in Ω B \Ω BA , any such reactions are implemented according to the chosen microscopic reaction method [43][44][45].

Simulation specifics
The Gillespie SSA [17] is used to simulate the above-described reactions in Ω BA , as well as the diffusive fluxes over the interface. The SSA requires the computation of an exponential random variable which gives the time, τ , until the next event, and can be found by transforming a uniform random variable u ∼ Unif(0, 1) by the following equation: Here, α 0 (t) is the sum of all of the propensity functions: and The PDE solutions and Brownian dynamics are implemented using the same discrete time-step, t, and the diffusive jumps across the interface (and any required reactions, r ∈ R * ) are implemented in an event-driven manner, according to the Gillespie SSA. Event-driven time-steps are implemented until the putative time for the next event passes the next Brownian/PDE update time, at which point the PDE and Brownian dynamics are updated. Pseudocode for the ARM is given in algorithm 1. , and their sums according to equations (3.14) and (3.15). Calculate α r (t), for r ∈ R * , using equation (3.11) and finally compute α 0 (t) according to equation (3.13). (1c) Calculate the time, τ , until the next auxiliary region event according to equation (3.12).
(ii) If u 1 α 0 (t) < α 0 P (t) (corresponding to a jump from Ω PA to Ω BA ): • Use u 2 to determine the species, k, which the jump affects, with each species selected with probability proportional to its propensity function. • Remove a particle from the PDE auxiliary region for species k via equation (3.9). • Initialize a new particle of species k within Ω BA at position y * = u 3 h a + I.
Else if u 1 α 0 (t) < α 0 P (t) + α 0 B (t) (corresponding to a jump from Ω BA to Ω PA ): • Use u 2 to determine the species, k, which the jump affects, with each species selected with probability proportional to its propensity function. • Choose a particle of species k uniformly at random from within the Brownian auxiliary region and remove it from the system. We do this by selecting an index q such that q = u 3 N k BA , where x denotes the smallest integer larger than x. • Add a new particle into the PDE auxiliary region for species k via equation (3.8).
Else (corresponding to a reaction in Ω BA ): • Use u 2 to choose the reaction r ∈ R * (t) to be implemented with probability proportional to its propensity function. (iii) Implement any reactions using an appropriate method (see §3.1). Note that production reactions should be implemented after any degradation reactions in order to prevent particles being created and destroyed in the same time-step. (iv) Set t = t , update t = t + t.

Results
Within this section we present four test problems which are used to demonstrate that the ARM correctly simulates reaction-diffusion systems. Two of these problems are models of pure diffusion with different initial conditions and will demonstrate that the fluxes over the interface are consistent with the expected behaviour of the fully Brownian simulations. The third problem is the formation of a morphogen gradient, which demonstrates the successful implementation of reactions in the ARM. Despite the fact that our method is valid for higher-order reactions, the first three test problems consider reactions up to first order. For such systems, no moment closure assumptions are required to derive the mean-field reaction-diffusion PDE and hence its behaviour agrees with the mean behaviour of the individual-based models. This allows us to efficiently verify accuracy by comparing the mean behaviour of our hybrid method to the known mean-field behaviour. Finally, in test problem four, we implement a second-order reaction system in higher dimensions, indicating the applicability of the method to more complicated examples.
For each of the first three test problems, we use Ω P = (−1, 0) and Ω B = (0, 1), meaning that the interface is the single point at 0. We take the value of the fixed PDE and Brownian time update steps to be t = 0.02, the auxiliary regions have width h = 0.05 and the diffusion constant is D = 0.0025 (unless specified otherwise). We will quantify the qualitative comparisons, presented in this section through  density comparison snapshots, in §5. All simulations will comprise only a single species, so henceforth, all sub-or superscripts, k, pertaining to species will be removed.

Test problem 1: maintaining equilibrium
For the first test problem, we simulate pure diffusion in the form of a simple Brownian motion with reflecting boundary conditions, which has Fokker-Planck equation given by the diffusion PDE and corresponding boundary conditions: where p 0 (x) denotes the initial condition. Note that p(x, t) here represents the mean-field solution across the whole domain, whereas c(x, t) represents the PDE solution in Ω P in the hybrid method. We initialize particles uniformly across the computational domain, so that p 0 (x) ≡ N/2, where N is the (constant) number of particles in the system. Figure 5 shows that the ARM passes the most basic test by maintaining the steady state without causing an accumulation of mass on either side of the interface. For this test problem, we also include a plot which displays the variance in the density of particles (figure 6). In order to calculate this variance, we have binned the positions of particles on the spatial domain onto a mesh of size h a (the same as the auxiliary region width) and calculated the variance of the density in each bin over a number of identically initialized (up to random allocation of particles in Ω B ) repeats. This demonstrates a problem that occurs with all hybrid methods which contain an interface coupling a stochastic to a deterministic region. The variance is damped close to the interface in the stochastic part of the domain, due to the deterministic nature of the solver on the opposite side. Specifically, the PDE effectively has a stochastic boundary condition at the interface, caused by the diffusive jumps between the auxiliary regions. This causes a higher level of variance than would be expected if it was a purely deterministic regime. However, when a particle jumps from the PDE to the Brownian dynamics auxiliary region, because the PDE region is mostly deterministic, it contributes less variance than would be expected than if the stochastic method was employed across the entire domain. There are methods that can be used in order to fix this problem, such as the use of an overlap region (e.g. [36]) and replacing the PDE with an appropriate stochastic PDE (SPDE) (e.g. [42]). This is explored in more detail in the discussion ( §6).     all particles uniformly within the PDE domain, Ω P , which results in

Test problem 2: flux over the interface
The results from this simulation are displayed in figure 7. As with the uniform initial condition in test problem 1, we see from figure 7 that the hybrid method agrees with the solution of the mean-field model, indicating that the method simulates flux over the interface accurately. We have also tested our hybrid method with all the mass initialized uniformly across Ω B and found a similarly good agreement between the hybrid method and the mean-field solution (figures not shown).

Test problem 3: morphogen gradient
For the third test problem, we investigate the formation of a morphogen gradient from a uniform initial condition. The gradient is formed by allowing particles to diffuse throughout the domain as well as to degrade at a rate μ. We also have particles entering at the left-hand boundary, For the corresponding microscopic dynamics, we implement Brownian motion for the diffusion of particles and a time-based method in order to enact the degradation reactions. We note that production of particles is not implemented within the microscopic domain as it occurs at x = −1. N = 500 particles are initialized uniformly across the whole domain [−1, 1]. As demonstrated in figure 8, the solution of the hybrid method matches that of the corresponding mean-field model, as with the previous two test problems.

Test problem 4: higher-order systems
For our final test problem, we look at the reaction system (4.5) which takes place in a three-dimensional cuboid Ω ⊆ R 3 of volume V, where Ω = (x 0 , x 1 ) × (y 0 , y 1 ) × (z 0 , z 1 ). We further divide this domain by firstly defining the position of the adaptive planar interface I(t) ∈ (x 0 , x 1 ) which is to be employed in this test problem (see §4.4.3). In an analogous way to in the onedimensional case, we then define the now time-dependent PDE and individual-based subdomains, Ω P (t) and Ω B (t), with volumes V P (t) and V B (t), respectively. These subdomains and volumes depend on t due to the adaptive interface position. The interface will move according to the local density profile within the PDE and Brownian dynamics auxiliary regions Ω PA (t) and Ω BA (t), which are explicitly defined to be and Ω BA (t) = (I(t), I(t) + h a ) × (y 0 , y 1 ) × (z 0 , z 1 ).
Before specifying how the interface will move, we will first find a PDE in one dimension that we will use to simulate the deterministic part of our domain. We do this by considering the reaction system (4.5) and forming an ODE in three dimensions. We then include isotropic diffusion to obtain a three-dimensional PDE, and finally impose a constraint on the initial condition to simplify this to a one-dimensional PDE. We then briefly describe the process we use to evolve the individual-level behaviour, before introducing an adaptive interface. We will finish this subsection with the results of some simulations of this system. Note that from now on, we will drop the dependence on t for any of the subdomains, their volumes and the interface position for brevity, unless they are explicitly needed.

PDE model
We will use the CME for the reaction system in order to derive a PDE that approximates the system If we now define the kth central moment A k := ∞ n=0 n k p n , we can multiply the CME by n and sum over all n ∈ N 0 to yield the mean equation: The ODE (4.6) is currently exact, but depends on the second moment of A. Furthermore, the ODE for every moment of A depends on higher moments still-the system is not closed. In order to close the system, we follow Erban & Chapman [43] and apply Poisson moment closure, which implies the following: Applying the moment closure (4.7) to the ODE (4.6), and setting c = A /V P gives us the closed ODE Finally, including isotropic diffusion through the usual Laplace operator yields the three-dimensional PDE: We will enforce an initial condition which is translationally invariant in both the y and z coordinates, which means that the dynamics will remain translationally invariant for all time. As such, c is simply a function of x and t, and the dynamics can be represented by a one-dimensional equivalent of this PDE by implementing zero-flux conditions on all boundaries and using the following transformation: where L y = y 1 − y 0 and similarly for L z = z 1 − z 0 . This gives the following:

Individual-based formulation
We now turn our attention to the individual-based system. In order to simulate the three-dimensional individual-based model, we will follow the λ-ρ method [43]. In the context of this system, whenever two particles are within the reaction radius ρ, they react with a probability P λ , which is a function of the kinetic rate κ 1 , the time-step t, and the diffusion coefficient D. For more information on how P λ is chosen, we refer the reader to [43]. The zeroth-order reaction is completed by initializing a particle uniformly throughout the individual-based domain Ω B with probability κ 2 tV B , which we ensure is below 1 by choosing t to be sufficiently small.

Adaptive interface
Test problems 1-3 have been simulated using a static interface. However, this requires a priori knowledge of where the interface should be for all time. When the finer scale modelling regime is required in order to resolve a fixed area of space in more detail (for example, the region around ion channels [12]), the interface position may be known. However, if the purpose of the interface is to split regions of space in which there are high and low particle numbers in situations in which particle numbers change dynamically, a different approach is required. In this case, the interface (or interfaces) need to move with the density of particles to maintain the computational savings they are designed to provide. We now describe a method, adapted from Robinson et al. [29], which allows the interface to move adaptively. The interface at time t, which we shall denote by I(t), moves according to local particle numbers in the auxiliary regions around it. We set two thresholds β u > β l , and move the interface towards the PDE  subdomain if N PA (t) < β l and towards the individual-based subdomain if N BA (t) > β u (borrowing the notation from §3). The two threshold values are designed to prevent the interface from rapidly oscillating between two values, which is a possibility when β u = β l due to the stochastic nature of the system. We enforce that when the interface moves, it moves a distance h a , the width of the auxiliary region, in the chosen direction. If the interface moves towards the PDE subdomain (i.e. N PA (t) < β l ), we convert the PDE auxiliary region into particles, initializing each one uniformly across the new Brownian auxiliary region Ω BA . As N PA (t) is not necessarily an integer, we treat the fractional part (N PA (t) mod 1) to be the probability of initializing one extra particle within the newly formed individual-based region. We then scale the rest of the PDE subdomain to ensure that we conserve mass. During an interface movement towards the individual-based subdomain (i.e. N BA (t) > β u ), the Brownian auxiliary region is converted to PDE mass by initializing a density of N BA (t)/h a uniformly across the new PDE auxiliary region, Ω PA , created by moving the interface. For a more detailed description of a similar method, we direct the interested reader to Robinson et al. [29].

Results
We consider N particles initialized across Ω with a constant negative gradient so that the density of particles at position x 1 is equal to zero. This ensures that the interface will move as the dynamics progress. The results can be seen in figure 9, in which the hybrid method has been averaged over S = 1000 repeats. The hybrid density in the case of the moving interface is represented as yellow bars throughout the domain. This is because the interface position changes with each repeat, and so very few regions of space are solely represented by one or the other modelling paradigm over all repeats.
We can see good agreement between the hybrid method and the fully individual-based method throughout the domain, with the only discrepancy close to the left hand boundary at 0 caused by the difference between the PDE-and individual-based methods due to moment closure. We compare our hybrid method to the fully individual-based method here, in contrast to the PDE solution used in test problems 1-3, due to the inaccuracy introduced in the PDE by the moment closure required for the second-order reaction.

Error analysis
We have seen in §4 that the solutions provided by the hybrid method visually match the mean-field solution. Within this section, we quantify the difference between the solutions of these test problems. We compare the mass in the PDE and Brownian regions of the domain using the two methods. Separately we   compare the density profile across the whole domain using the histogram distance error (HDE). We then proceed to investigate the dependence of the accuracy of the hybrid method on the two free simulation parameters ( t and h a ).

Quantitative comparisons
In order to evaluate the accuracy of the ARM for test problems 1, 2 and 3, we compare its mean behaviour (averaged over S = 1000 repeat simulations) to the mean-field model for which we compute the analytical solution across the entire domain Ω, for each of our test problems. Figure 10 contains nine plots which demonstrate the error for the first three test problems above; figure 10a-c is for test problem 1, figure 10d-f is for test problem 2 and figure 10g-i for test problem 3. The first and second columns show particle number comparisons between the hybrid and analytical solutions. Specifically, in the first column we compare The last column of figure 10 contains the HDE, which is defined by This ensures a value of the HDE between 0 and 1. Here, 0 means that the two solutions are exactly the same, and 1 corresponds to the two solutions having non-overlapping supports. All figures were produced using the same number of repeats (S = 1000). In all cases, the relative errors between the mean-field and hybrid methods, in figure 10, are low with no discernible bias about zero. Similarly, all HDE plots in figure 10 are low for the majority of the simulations. This demonstrates numerically that the hybrid scheme presented in this paper is correctly reproducing the behaviour of the Brownian model in the mean field. These error plots confirm the visual concurrence shown in figures 5-8. For the fourth test problem, we use a different error measurement due to the disparity between the mean-field PDE and individual-based systems. Consequently, we choose to compare the number of particles in the final compartment for both the hybrid method and the individual-based method (figure 11). We motivate this in two ways. Firstly, using this measure of error, we are able to minimize the influence of the extra error caused by the difference between the mean-field PDE and the individual-based method. Secondly, several biological systems require detailed knowledge of the particle concentrations at the end of the domain. Apical growth of filamentous cells such as fungi [46] is such an example. If we define N H (t) to be the average number of particles in the region (x 1 − h a , x 1 ) × (y 0 , y 1 ) × (z 0 , z 1 ) (where we recall that Ω = (x 0 , x 1 ) × (y 0 , y 1 ) × (z 0 , z 1 )), when simulating the hybrid method at time t, and the quantity N M (t) to be the same for the fully microscopic simulation,  Figure 11. The error measurement given in equation (5.6) for the system simulated in figure 9.
we can obtain a measurement of error given by The relative error shows no long-term bias in either direction, and oscillates around zero, indicating a close agreement between our hybrid method and the ground truth individual-based method. The hybrid method completed 1000 repeats in 485.5 s, while the fully individual-based method took 1047.4 s.

Parameter choice
Within the ARM, there are two free parameters that need to be chosen-the width of the auxiliary regions h a and the time-step for the PDE and Brownian updates t. These need to be chosen so that the quantity D t/h 2 a remains small enough that the particle numbers in the auxiliary regions do not become overly equilibrated between PDE/Brownian update steps. That is to say, if there is a gradient across the interface, t should be small enough that the closed system of the two auxiliary regions should not reach steady state between PDE/Brownian update steps.
In order to demonstrate why D t/h 2 a must be small, we consider the evolution of particle numbers in the two auxiliary regions between PDE/Brownian update steps. We form an ODE for particle numbers in one of these boxes (using the fact that particle numbers are conserved between PDE/Brownian updates).
Let ν 0 be the (constant) number of particles in the two auxiliary regions combined, M P (t), M B (t) be the mean number of particles in the PDE and Brownian auxiliary regions, respectively, at time t, and μ P , μ B be the number in the PDE and Brownian auxiliary regions, respectively, at time 0, which will represent the beginning of a time-step. Then, the equation for the mean number of particles in the PDE auxiliary region can be calculated from a simple probability master equation as where we recall that d is the jumping rate between the two auxiliary regions and is linked to the diffusion constant, D, via equation (3.3). Solving this ODE gives Assuming a small time-step, t, we can approximate M P ( t), the number of particles after a time-step has occurred, by Taylor expanding equation (5.7) to first order:  Fixing the value of D and using equation (3.3), we find that We require the change in the number of particles over the small time-step to be small, and so would like M P ( t) ≈ μ P . Thus we need to choose our parameters such that the quantity D t/h 2 a small. This elucidates an important relationship between the fixed and free parameters of the model. If the diffusion coefficient is large, then we must choose a small update time-step or a larger auxiliary region length to compensate. Figure 12 shows that a large region of the t-h a space has a very low HDE, meaning that our method is robust to parameter change, and only breaks down once the value of D t/h 2 a becomes very large. The plot also shows that, given any choice of the width of the auxiliary regions, h a , there is a value for the time-step, t, which will give a low level of error. Also, depending on our choice of t, we can adjust h a to make the simulation more accurate.

Discussion
We have presented a new spatially coupled hybrid method for coupling a Brownian dynamics representation of a reaction-diffusion system to its corresponding mean-field PDE description. By bridging the gap in spatial scales with intermediate auxiliary regions, we have produced an algorithm that is not only accurate, but is also robust to the choice of the free parameters within the problem, namely the width of the auxiliary regions, h a , and the fixed time-step, t, used to update both the PDE and Brownian dynamics. This is in direct contrast to a previously presented PDE-to-Brownian hybrid, which we demonstrated to be extremely parameter-sensitive. In order to make the ARM even more robust, applicable and efficient, we now discuss several areas for possible extension, which will be addressed in future works.
In the interests of completeness we should point out that, as with the pseudo-compartment method of Yates & Flegg [37], the ARM requires that the mass in the PDE auxiliary region Ω PA be sufficient for a step function, corresponding to the mass of a particle, to be removed uniformly from across the auxiliary region. This will lead to difficulties in situations in which particle numbers are low around the interface. Arguably though, we should not employ such hybrid methods in situations for which particle density is low around the interface as the PDE will be a poor model of the true stochastic, microscopic dynamics in these regions. A possible solution to this inconvenience is the incorporation of an adaptive interface, which we have employed in test problem 4. Such interfaces evolve with the simulation dynamics, ensuring the appropriate model is used for the corresponding particle density [29].
A related issue is that of multiple interfaces. Multiple interfaces will allow the efficient simulation of stochastic reaction-diffusion systems in which multiple regions of high and low concentration are  expected. Such patterns will require interfaces to be dynamic in number and transient in nature.
Although we have not implemented such interfaces in this work, we expect it to be a relatively straightforward extension. While we have presented an example in which the system is simulated in a cuboid with a planar interface (test problem 4), non-planar interfaces, such as those which have corners or are curved, and complex domain geometries, present deeper challenges that we hope to address in a future publication. Failing to maintain stochastic variation is a problem which is common among many spatially coupled hybrid methods. As a result of the deterministic nature of the PDE, the noise in the Brownian dynamics region of the domain is damped in comparison to the fully microscopic model ( figure 6). In the literature, two approaches have been used in order to rectify this. The first is an overlap region, which has been employed in several papers [36,41,47]. These methods introduce a region of space which lies in the intersection of the two domains. In these regions, mass is simultaneously represented using both scales of description. The second is to replace the deterministic PDE with an appropriately chosen SPDE. Alexander et al. [42] consider such a coupling and demonstrate they can indeed fix the discrepancy by using an SPDE as their continuum macroscale model. We will address the use of both SPDEs and overlap regions (in which the region between the PDE and the Brownian dynamics regions is simulated using a purely compartment-based method) in forthcoming work.
The ARM provides a simple yet accurate method to couple an individual Brownian dynamics representation of a reaction-diffusion system to a corresponding PDE representation. Our hybrid algorithm will be of particular interest to researchers modelling reaction-diffusion systems whose concentrations vary significantly across the spatial domain. By reducing the computational expense of simulations, the ARM will facilitate the investigation of stochastic effects in such systems, in some cases, making the difference between being able to interrogate the system and not. In particular, we suggest that our method will be useful for the investigation of stochastic Turing patterns [48], Fisher waves [30,31], oscillatory dynamics [49] and excitatory dynamics [50] with applications in embryogenesis [4], intracellular dynamics [6] and pattern formation [48] among others. It may also be worthwhile to interface the methods presented here with commonly used Brownian dynamics simulation software packages such as Smoldyn [9].   Figure 14. The error plots for the example in figure 13. The first two panels show relative errors in particle numbers in the (a) PDE and (b) Brownian dynamics subdomains. The histogram distance error is displayed in (c).
parameter that is to be defined is the auxiliary region width, which we set here to be h a = 0.1. The results can be seen in figure 13.
As can be seen from this figure, the agreement between the mean-field and hybrid solutions is much closer than that of the PDE-assisted Brownian dynamics [41]. This indicates an improvement over the previous method. We also present the error plots which are described in §5-namely the relative errors in particle numbers and the HDEs.
Once again, the relative error plots (figure 14a,b) appear to show no long-term bias in either direction and the HDE (figure 14c) is small.