{'fraction_non_alphanumeric': 0.025762980578676178, 'fraction_numerical': 0.006089431773141642, 'mean_word_length': 4.68846075015372, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A hybrid algorithm for coupling PDE and compartment-based dynamics

Stochastic simulation methods can be applied successfully to model exact spatio-temporally resolved reaction-diffusion systems. However, in many cases, these methods can quickly become extremely computationally intensive with increasing particle numbers. An alternative description of many of these systems can be derived in the diffusive limit as a deterministic, continuum system of partial differential equations. Although the numerical solution of such partial differential equations is, in general, much more efficient than the full stochastic simulation, the deterministic continuum description is generally not valid when copy numbers are low and stochastic effects dominate. Therefore, to take advantage of the benefits of both of these types of models, each of which may be appropriate in different parts of a spatial domain, we have developed an algorithm that can be used to couple these two types of model together. This hybrid coupling algorithm uses an overlap region between the two modelling regimes. By coupling fluxes at one end of the interface and using a concentration-matching condition at the other end, we ensure that mass is appropriately transferred between PDE- and compartment-based regimes. Our methodology gives notable reductions in simulation time in comparison with using a fully stochastic model, whilst maintaining the important stochastic features of the system and providing detail in appropriate areas of the domain. We test our hybrid methodology robustly by applying it to several biologically motivated problems including diffusion and morphogen gradient formation. Our analysis shows that the resulting error is small, unbiased and does not grow over time.

Multiscale modelling challenges occur frequently throughout cellular biology and in the context of cell migration. Spatial reaction-diffusion models can be used to describe, either deterministically or stochastically, various biological phenomena. These include actin dynamics in filopodia 1 , calcium signalling 2 and chemisorption of polymers 3 . In many cases, it may be beneficial to use a multiscale approach to modelling using different descriptions in different spatial regions. In this article, we will set out a method for coupling together a continuum deterministic description and a discrete stochastic description.
Commonly, continuum approaches using partial differential equations (PDEs) are adopted to model biological systems 4 . These equations can either be solved analytically (in some cases) or simulated numerically. Results using this methodology are relatively fast to calculate computationally. However, for systems with small numbers of molecules the results obtained using deterministic methods may not always capture the behaviour of a stochastic system appropriately, especially in situations where molecular numbers are low and interactions are non-linear. For example, in a system with multiple steady states such as the canonical model of Schlögl 5 , a deterministic model fails to capture the switching behaviour between the steady states seen in a stochastic model. In general, PDE models break down when the number of molecules present is very low and stochastic effects dominate 6 . Although deterministic models may provide useful information about average behaviour (in the case of linear systems) they cannot offer a full description of every system. Thus, in the case where copy numbers are low, the best description will be afforded by a stochastic model. There are two main types of stochastic models used for reaction-diffusion equations 7 : off-lattice methods and on-lattice compartment-based methods. We will focus on compartment-based methods, which generally offer a coarser description than their off-lattice counterparts. When simulating a system using a compartment-based stochastic model (also known as a position-jump model), the computational cost of the simulations can become prohibitive if the number of particles in the system is high. A computationally efficient continuum model may be more appropriate in this scenario. Thus in situations where particle concentrations vary widely across the domain there may be advantages to using a continuum PDE model in the region of the spatial domain where particle numbers are high and a discrete stochastic model elsewhere. Moreover, detail is often only required in a certain part of the domain and thus a spatial-hybrid model may be most appropriate 1,3,[7][8][9][10] . Such a hybrid model would allow an accurate representation of the reaction-diffusion dynamics in the region where this is required but minimises the computational resources needed to perform the calculation by using less detailed, more efficient methods in regions of the domain where detail is not required.
Previously, Flekkøy et al. 8 have developed a hybrid model that links a PDE-based model to the motion of random walkers on a lattice. Motivated by heat transport around a facture in a solid, Flekkøy et al. 8 choose a detailed description of the particle dynamics coupled to a coarse-grained PDE model: the lattice spacing used to solve the PDE is larger than in the particle-based region. More recently, PDE-to-compartment hybrid methods have been developed which employ a region of the PDE regime in which particles are represented using both the compartment-and PDE-based modelling regimes simultaneously 11,12 . The duality of these so-called "pseudo-compartment methods" allows for particles to behave correctly as they cross individually between the two different regimes since particles can jump into their neighbouring compartment according to standard compartment-based rules for diffusion.
Our hybrid modelling regime employs a PDE mesh that is significantly (and arbitrarily) finer than the lattice in the compartment-based region. This choice is natural in many situations, including in a biological context, where we are choosing to use the PDE model in regions of high population to offer improved computational efficiency. Taking a fine mesh will not prove computationally prohibitive compared to the stochastic model, but allows us to make the numerical solution of the PDE arbitrarily accurate. Methodologies with coarser or equal PDE spacing relative to compartment spacing 8,12 are open to questions about what exactly the "PDE regime" represents given its resolution and accuracy are restricted by the resolution of the compartment-based method.
Our approach to coupling of deterministic PDE-based and stochastic compartment-based regions employs an overlap region where both modelling descriptions are valid. This overlap region can contain multiple compartments if desired. The method that we have developed relies upon specifying a Dirichlet-type condition between the two models at one interface at the edge of the PDE-based region and dictating the correct flux of PDE on compartments at the other interface. This fixes the boundary conditions at the interfaces between each of the regions.
In the remainder of this article, we describe and explore our novel hybrid coupling al-gorithm in detail and illustrate the effectiveness of the method. In section I B, we present the hybrid method in full and justify the coupling conditions chosen. Thereafter, in section II A, we demonstrate the appropriate behaviour of our method through its application to systems of diffusing particles with various extreme initial conditions (chosen specifically to test the algorithm) and a biologically-motivated example: the formation of a morphogen gradient. We apply the model to a travelling wave example in section II A 3 and introduce an adaptive interface between the modelling regimes in section II A 4. We present detailed simulation-time comparisons of the hybrid model with the fully stochastic model for our test problems at the end of section II A 4, which explicitly demonstrate the improved efficiency of our hybrid method. The fidelity of the algorithm's performance is then examined and the error (with respect to a range of model parameters) analysed in section II B. We conclude in section III, with a discussion of the potential advantages of this hybrid method in relation to other existing methods.

A. The domain
Suppose, arbitrarily, we have a domain Ω = [−1, 1] which we divide into a region Ω c in which we use a compartment-based, stochastic model and a region Ω p in which we use a deterministic, PDE-based model. A characterising feature of our hybrid methodology is an overlap region (shown in Figure 1) in which both modelling regimes are simultaneously valid descriptions (i.e. Ω c ∩ Ω p = ∅). Either side of the overlap region, we have interfaces I 0 and I 1 (see Figure 1). In a similar context, it has been demonstrated that an overlap region is required to give the appropriate variance for a model coupling a Brownian motion particle-based description and a PDE-based model 3 .
In Ω c , the domain is split into compartments of width h, where the kth compartment  sufficiently high particle numbers. For low concentrations, this assumption breaks down. If the concentration of particles is u(x, t), given a total of N particles in the system, then we can relate the probability of finding any particular particle in the system, p p (x, t), to the concentration as p p (x, t) = u(x, t)/N. This probability density remains well defined even at low particle copy numbers, when we cannot interpret the PDE as a concentration but are able to view it as a probability density. Therefore we interpret p p (x, t) as the probability density to find each of N particles in the PDE region at a given location x at a certain time t. The expected number of particles in a subset, ω, of the PDE domain, Ω p , is given by N ω p p (x, t)dx. We will use u k (t) where k = 1, ..., K + 1 to denote the PDE density at the kth PDE lattice point in the finite difference discretisation of the PDE required for our hybrid algorithm.
For the compartment-based regime, let p c (x, t) (defined initially only at the centre of compartments) be the probability of finding one of the identically initialised particles at position x at time t. Since each compartment is well mixed, we can describe the evolution of p c (x, t) using the reaction-diffusion master equation 14 . We will also use the notation .., A n (t)) to represent the distribution of particle numbers across compartments.

B. The coupling algorithm
We now describe an algorithm which couples the two regimes together. Informally, the coupling is achieved by setting the value of the PDE lattice point at I 0 to the average of the adjacent compartment populations in Ω c and using the gradient in the PDE-based region Ω p to give a rate of jumping across the interface I 1 for the compartment-based regime.
In what follows, we specify and justify these coupling conditions mathematically. These conditions are analogous a to Neumann condition for the compartments at I 1 and a Dirichlet condition for the PDE at I 0 . We aim to apply an appropriate flux of particles to and from Ω c based on the PDE profile across the interface I 1 , which will ensure that the gradients of the different modelling regimes agree. Feasibly, if this were the only condition, situations could arise where the gradients of the two regimes agree, but there is a notable discontinuity in the values of the density between descriptions. To prevent this, we enforce a boundary condition on the PDE requiring the density on the lattice point at I 0 to match an average of the density of the surrounding compartments.
First consider the Dirichlet matching condition at I 0 , where we specify the PDE density in terms of particle numbers: Writing this in terms of the analogous probability densities, we have: Extending p c to continuous space and Taylor expanding the terms on the right-hand side (RHS) of equation (2) to first order, we find that: This suggests that matching condition (1) ensures agreement between the solution of the PDE and compartment-based particle numbers at I 0 . The agreement will become exact in the limit h → 0.
For the condition at I 1 , we want to match the flux across I 1 in the compartment-based regime to that in the PDE regime. We will show that by enforcing the matched flux condition, the probability density for the compartment-based region evolves according to the diffusion equation in the limit of small compartment size.
We begin by writing down the mater equation 14 for the probability density of a single particle at compartment −m, adjacent to the interface, I 1 : where p k c (t) is shorthand for p c (I 0 + (2k + 1)h/2, t) and describes the probability density for a single diffusing particle to be found in the kth compartment at time t. Here, ψ p is the flux imposed (as part of the hybrid algorithm) on compartment −m from the right. If there were compartments to the right of the compartment labelled −m (i.e. −(m + 1) etc (see Instead we must approximate the true flux, ψ c , by an ansatz derived from the PDE, ψ p as follows.
Suppose that the lth lattice point of the PDE lies on the interface I 1 , and w is the ratio of spacing between the compartment size, h, in Ω c and the PDE finite difference lattice size, We can interpolate the density in Ω p at the centre of the −mth compartment by where Imagine an extra compartment −(m + 1) to the right of I 1 . We could interpolate the density at the centre of this ghost compartment using a similar expression: where Given these interpolations of the PDE density at the centre of compartments we can approximate the diffusive flux across the interface and consequently set Substituting this into equation (3) gives Upon rearrangement this implies In order to demonstrate the veracity of our choice of ψ p , we extend p c to be a continuous function of space and Taylor expand terms on the RHS in space about the centre of the −mth compartment (i.e. I 1 − h/2). Taylor expanding p c I 1 − h 2 , t + ∆t in time and taking the diffusive limit, we find we recapitulate the diffusion equation for the probability density Consequently, this indicates that the flux ψ p given by equation (5) is an appropriate boundary condition for the compartment-based model. Given our two matching conditions at either end of the interface, the hybrid algorithm can be implemented in a time-driven sense as follows: (i) Initialise number of particles A k , k = −m, ..., −1, 1, ..., n in compartments in Ω c and apply consistent initial conditions in Ω p .
(ii) Select the compartment-based time step ∆t, such that the probability of more than one event occurring per time-step is O(∆t 2 ), and a maximum duration of the simulation, T f inal . Set t := 0.
Draw a uniform random number r 1 .
(iv) Calculate a uniform random number r 2 .
If r 2 < α 0 ∆t where α 0 is the total propensity of the 'reaction' events in the compartment-based regime, then a reaction occurs in that time step.
(v) If a 'reaction' occurs, generate a uniform random number r 3 , and find j such that Update number of particles in each compartment according to chosen reaction, j.
(vi) Update time such that t := t + ∆t.
(vii) Update PDE region Ω p using an appropriate numerical method. Apply the boundary condition at the right-hand boundary and the coupling condition at I 0 as follows: (viii) If t < T f inal , then go back to step (iii). Else end.
Algorithm 1: Time-based hybrid algorithm for stochastic reaction-diffusion simulations using a compartment-based region and an overlapping PDE-based region.
Note the factor of N in the calculation of ψ at step (iii) is due to the scaling between concentration and the probability distribution for a single particle. Both "time-based" and "event-based" versions of the hybrid coupling algorithm are possible. The main difference between these is that the time-based algorithm uses a fixed time step ∆t to update both Ω c and Ω p , while the event-based algorithm steps forward to the next reaction in Ω c , while still fixing a maximum time step in Ω p for updating the PDE. For systems with large numbers of particles, the event-based algorithm will be more efficient as it allows the use of larger time steps in the stochastic regime so fewer steps of the algorithm are required. However, for simplicity, we present here the time-based version, Algorithm 1.

Test problem: diffusion
We will begin our examination of practical applications of the hybrid coupling algorithm by applying the method to a test problem in which particles diffuse with diffusion constant D. With large copy numbers of particles in the system, the density of diffusing particles, u(x, t), is governed by the diffusion equation: Adding reactions to this system should not affect the boundary behaviour directly and therefore it is sufficient to test our model on a problem of this type 15 . As previously specified We consider three different initial conditions, φ(x): a uniform initial condition, demonstrating that the algorithm can maintain an equilibrium state, a step function with all the mass in [0, 1], that is: and a step function with all the mass in [−1, 0], that is: These provide a robust test of our hybrid algorithm in a variety of different scenarios, showing it can maintain net flux from each region to the other.
We have performed simulations of the hybrid model, using the three different initial conditions described above. We also present the analytical solutions of the mean-field diffusion equation. In particular, suppose that all the mass is initially in [0, 1], as in (8). Using a Green's function and an infinite series of images at the boundaries we obtain an analytical solution to equation (7) of the form: where we have written the solution in terms of error functions. The solution for the initial condition of a step function with all the mass in [−1, 0] (as in equation (9)) can be obtained by symmetry from equation (10). This solution is used in both Figure 3 and

Test problem: morphogen gradient
We also apply our model to another test problem: the formation of a morphogen gradient.
For this problem, we use the same domain and partitioning as before. Morphogen molecules are produced at rate J at x = 1 and throughout the domain morphogen molecules decay with constant rate µ and diffuse with diffusion coefficient D. When there are sufficiently many molecules in the system, we expect the density of molecules, u(x, t), to be governed by the following PDE:  We apply zero flux conditions at the boundaries and initially we assume there are no molecules in the system.
The results of simulating this morphogen system are shown in Figure 6. The system was simulated up until t = 20 after which point the system had approached steady state. Good agreement can be seen between the hybrid simulation algorithm and the analytical solution of (11).

Test problem: travelling wave
The occurrence of travelling waves is common throughout the natural world: they describe a variety of phenomena from propagation of genes in a population 16 , to epidemic outbreaks 17 , and in the FitzHugh-Nagumo equations for a nerve axon pulse 18 . One commonly used model for a travelling wave front is the Fisher-KPP equation: where D is the diffusion coefficient, and k 1 and k 2 are reaction rates. This is a non-linear reaction-diffusion equation for the concentration or population density u in one dimension.
It can be shown that this results in the formation of a travelling front with a minimum wave-speed of c = 2 √ Dk 1 , given continuous initial conditions with compact support 19 .
Consider the reversible chemical reaction: Using the law of mass action in a deterministic setting 19 and including diffusion effects, results in the Fisher-KPP equation (12) where v min is the minimum velocity for the deterministic model, K is a constant, and N is the total number of molecules 21 . By considering moments of the appropriate chemical master equation, we obtain a hierarchy of coupled equations, where the kth moment depends upon the (k + 1)th moment 22 . In order to obtain a closed system we must make a closure approximation. The degree of agreement between the deterministic and stochastic descriptions will depend on the validity of this closure assumption.
We note that the nature of reaction scheme (13) means that population growth in compartments ahead of the wave front does not begin until there is at least one particle present in that compartment. The discretisation of particles in the stochastic model, therefore, restricts the progress of the wave and results in the lower wave speed in comparison to the deterministic interpretation 23 .
Given that we do not expect the stochastic model to correspond to the deterministic model in the mean-field we will use a fully stochastic compartment-based description of the system for comparison with our hybrid system in order to determine its accuracy (as opposed to the PDE description which represented the mean-field behaviour of the previous test systems). We expect to make computational savings by using a PDE to describe the mean field behaviour behind the wave whist using the stochastic compartment-based model to simulate behaviour at the wave front and ahead of the wave, which determines the wave speed.
Applications of hybrid models to travelling waves have been made in previous work.
Moro 9 have successfully demonstrated such a model, using a flux based approach similar to that of Flekkøy et al. 8 . This hybrid model was then used to confirm the scaling of the velocity correction for the stochastic mesoscopic model. Further to this, an adaptive version of the two-regime method has also been applied to a travelling wave problem 20 . This model couples a microscopic Brownian motion based description to a mesoscopic compartmentbased description, as in the original two regime method 7 . In addition, the interface between the two regions is, in this case, allowed to move adaptively following the propagation of the front 20 . This enables the microscopic description to represent the most appropriate region of the domain, following the front of the wave, with the less computationally intensive mesoscopic description remaining behind the wave.
We demonstrate that our hybrid model can be applied successfully to a travelling wave An important measure when investigating stochastic simulations of reaction system (13) is the resulting wave speed. It can be difficult with a stochastic model to specify exactly where the wave front is at a given time and to quantify exactly how fast it is moving, since there will inevitable be noise in the results of simulations 22 . We choose to use the method outlined by Robinson et al. 20 , which considers the rate of change of the total mass, M(t) in the system. For times t 2 and t 1 , we take: where the factor k 2 k 1 is necessary since the height of the wave will approach k 1 k 2 . Dividing the rate of change of mass by this factor gives a measure of how fast the wave is propagating through the domain. between the wave-speeds of the two models as estimated by a moving average of the wavespeed estimates, after an initial transient. The slower initial wave-speed observed in both models is explained by the steep initial condition, which first needs to approach the profile of the travelling wave before it starts to move at constant speed.

Adaptive interface via a local detection criterion
In certain situations, as with the travelling wave presented in the previous section, the region of interest with lower particle numbers changes position dynamically. In order to capture most effectively the detail in this area whilst reducing the computational requirements it will be useful to have an interface that also changes position, so that regions with higher particle numbers can more often be modelled using the PDE. To ensure that the interface moves correctly, we initiate moves of the interface adaptively based on a detection condition of the particle density near the interface. Moving the interface should also prevent unnecessary simulation of large particle numbers using the stochastic regime in regions where this is not required; for example, behind the wave front for larger times in the travelling wave model (see Figure 7 panel (c)).
Such adaptive interfaces have previously been implemented in hybrid models 12  In the particular case of the travelling wave considered in the previous section (II A 3), it is important that we keep the entire front of the wave in Ω c , since it is the description governing the wave front that dictates the wave speed. To ensure this, we take u max = 10.5, δu = 1.0 for model parameters as in Figure 9.
When we have performed several iterations of the hybrid adaptive algorithm and wish to take an average of the results we encounter some difficulties. After a full iteration of the algorithm has been completed, the interface between the models will have, in general,   If I 0 has been updated, then density in newly created region is equal to density of that region in previous description.

Return to step (ii) unless final time is reached.
Algorithm 2: Algorithm for stochastic reaction-diffusion simulations with an adaptive interface using a compartment-based region and a PDE-based region.
Notable computational improvements are afforded by the hybrid model in comparison to the fully stochastic compartment-based model. Simulation time is decreased by a simulationdependent factor of around 5. Note that the adaptive interface algorithm for the travelling wave simulations is significantly faster than the scenario with the fixed interface.

Model
Fully

B. Sensitivity analysis
We demonstrate robustness of the coupling algorithm to choices of the algorithm parameters h, the compartment width, and ∆x p the PDE discretisation, showing how the total error varies as a function of these parameters. Since we are also able to vary the size of the overlap region in our coupling algorithm, we also demonstrate the effects of varying the number of compartments in this region. As the test problem here, we use simple diffusion with the same step-function initial condition as in Figure 3 given by equation (8). The results are presented in Figure 10.  to cell biology such as calcium signalling 2 , and wound healing 25 . In particular we focused on systems with multiple scales where detailed modelling is required in a certain region, but it might prove computationally wasteful to apply that method throughout the domain. Such systems occur frequently in a biological context due to the multiscale nature of biological systems 26 .
The hybrid algorithm that we have developed was robustly tested and demonstrated by applying it to several biologically motivated problems in section II A. There are noteworthy improvements in simulation time in comparison to a fully stochastic model, including a decrease in simulation time by approximately a factor of 5 when applied to a suite of standard test problems. The performance of this hybrid algorithm and the error compared to an analytic solution was analysed and explained.
At low particle numbers, a deterministic modelling method may no longer be appropriate and a stochastic method should be applied to account for the variation. There are disadvantages to the stochastic methods too; in particular they can require long simulation times.
In order to make best use of the complementary advantages of deterministic and stochastic models, multiscale hybrid models are becoming increasingly widespread, particularly in applications relating to reaction-diffusion systems. We have presented our own hybrid cou-pling algorithm to segue between stochastic compartment-based models and deterministic