Multiscale Stochastic Reaction–Diffusion Algorithms Combining Markov Chain Models with Stochastic Partial Differential Equations

Two multiscale algorithms for stochastic simulations of reaction–diffusion processes are analysed. They are applicable to systems which include regions with significantly different concentrations of molecules. In both methods, a domain of interest is divided into two subsets where continuous-time Markov chain models and stochastic partial differential equations (SPDEs) are used, respectively. In the first algorithm, Markov chain (compartment-based) models are coupled with reaction–diffusion SPDEs by considering a pseudo-compartment (also called an overlap or handshaking region) in the SPDE part of the computational domain right next to the interface. In the second algorithm, no overlap region is used. Further extensions of both schemes are presented, including the case of an adaptively chosen boundary between different modelling approaches.


Introduction
Stochastic models of well-mixed chemical systems are traditionally formulated in terms of continuous-time Markov chains, which can be simulated using the Gillespie stochastic simulation algorithm (SSA) (Gillespie 1977) or its equivalent formulations (Gibson and Bruck 2000;Cao et al. 2004;Klingbeil et al. 2011). These algorithms provide statistically exact sample paths of stochastic chemical models described by the corresponding chemical master equation (CME). However, they can be computationally expensive for larger chemical systems, because they explicitly simulate each occurrence of each chemical reaction. A number of approaches have been developed in the literature to decrease the computational intensity of SSAs. Taking into account separation of timescales, chemical reaction networks can be simplified by model reduction before they are simulated (Kang 2012;Kang and Kurtz 2013;Kang et al. 2014;Kang et al. 2019). The idea of model reduction can also be used to develop computational methods which efficiently estimate quantities of interest from stochastic simulations (Cao et al. 2005a, b;Erban et al. 2006;Cotter et al. 2011). Another approach is to describe the molecular populations in terms of their concentrations that change continuously (rather than treating them as discrete random variables). This can be achieved by the chemical Langevin equation, which is a stochastic differential equation (SDE) acting as a bridge between discrete SSAs and deterministic reaction rate equations (Kurtz 1976(Kurtz , 1978Gillespie 2000). Efficient algorithms which make use of the SDE approximations have been developed for the simulation of chemical systems especially when they include processes occurring on different timescales (Haseltine and Rawlings 2002;Salis and Kaznessis 2005;Griffith et al. 2006;Cotter and Erban 2016). More recently, the SDE approximations have been extensively used to develop hybrid algorithms which use boths SSAs and SDEs for different components of the studied systems (Liu et al. 2012;Ganguly et al. 2015;Duncan et al. 2016;Altintan et al. 2016). The chemical Fokker-Planck equation corresponding to the chemical Langevin equation can also be used to efficiently estimate quantities of interest from stochastic models Cotter et al. 2013;Liao et al. 2015;Cucuringu and Erban 2017).
In this paper, we consider spatially distributed (reaction-diffusion) models which can be described in terms of the reaction-diffusion master equation (RDME) (Erban et al. 2007). A spatial domain is discretized into compartments (which are assumed to be well mixed), and diffusion is modelled as a jump process between neighbouring compartments Kang et al. 2012b;Hu et al. 2014). In the literature, the RDME approach has been adapted to model and simulate spatially distributed systems using uniform meshes (equivalently, subvolumes or compartments) (Stundzia and Lumsden 1996;Elf et al. 2003;Lampoudi et al. 2009), nonuniform meshes (Bernstein 2005) or complex geometries (Isaacson and Peskin 2006). The resulting compartment-based model can be simulated by the Gillespie SSA. Compartment-based reaction-diffusion approaches have been used to model several intracellular processes, including Min oscillations in E. coli (Fange and Elf 2006;Arjunan and Tomita 2010), ribosome biogenesis (Earnest et al. 2015(Earnest et al. , 2016, actin dynamics in filopodia (Zhuravlev and Papoian 2009;Erban et al. 2014) and pattern formation in morphogen signalling pathways (Kang et al. 2012a). They have also been implemented in a number of software packages including MesoRD (Hattne et al. 2005), URDME (Engblom et al. 2009), STEPS (Wils and De Schutter 2009), SmartCell (Ander et al. 2004), Lattice Microbes (Roberts et al. 2013) and Smoldyn (Robinson et al. 2015). As in the case of the simulation of well-mixed systems, the Langevin approach provides an approximation of the compartment-based model which can reduce the computational intensity of simulations. Spatial Langevin approaches (Kalantzis 2009;Ghosh et al. 2015;Bhattacharjee et al. 2015) and stochastic partial differential equations (SPDEs) (Dogan and Allen 2011;Alexander et al. 2002Alexander et al. , 2005Atzberger 2010;  have been suggested to model stochastic reactiondiffusion systems. A hybrid method has also been introduced using the Langevin approximation for diffusion coupled with the compartment-based model for reactions (Lo et al. 2016).
In the thermodynamic limit (of large populations), compartment-based models lead to reaction-diffusion partial differential equations (PDEs) which are written in terms of spatio-temporal concentrations of chemical species. This property can be exploited to design multiscale (hybrid) algorithms which use the compartment-based Markov chain model in a subset of the simulated system and apply reaction-diffusion PDEs in other parts (Kalantzis 2009;Ferm et al. 2010;Yates and Flegg 2015;Spill et al. 2015;Harrison and Yates 2016). Other hybrid methods have also been developed in the literature including methods which couple more detailed Brownian dynamics (molecular-based) approaches with the compartment-based method (Flegg et al. 2012;Klann et al. 2012;Flegg et al. 2015;Dobramysl et al. 2016) or with reaction-diffusion PDEs (Franz et al. 2013;Schaff et al. 2016;Bakarji and Tartakovsky 2017).
In this paper, we analyse two multiscale algorithms which couple compartmentbased models with suitably discretized SPDEs. They can be used when a large number of molecules of some species are located in parts of the computational domain. In the region with a small number of molecules, we use a compartment-based model written as a continuous-time Markov chain. In other regions, we use SPDEs derived from the Markov process. The goal of this multiscale methodology is to get an approximation of the spatio-temporal statistics which we would obtain by running the underlying Markov chain model in the entire computational domain. The paper is organized as follows. In Sect. 2, we present the derivation of the SPDE description from the compartment-based model. In Sect. 3, two multiscale schemes are presented. An illustrative example with a static boundary between the SPDE and Markov chain subdomains is studied in Sect. 4. The algorithm is extended to a time-dependent interface in Sect. 5. In Sect. 6, we discuss an example with multiple species.

From Continuous-Time Markov Jump Processes to Stochastic Partial Differential Equations
We consider a system of N chemically reacting species S 1 , S 2 , …, S N , which are diffusing (with diffusion constants D i , i = 1, 2, . . . , N ) in the bounded domain Ω ⊂ R 3 . We use a compartment-based stochastic reaction-diffusion model (Erban et al. 2007), i.e. we divide the domain Ω into K compartments C k , k = 1, 2, . . . , K , and model the diffusion as a jump process between neighbouring compartments. In order to simplify the analysis, we consider that Ω is an elongated pseudo-one- Fig. 1a. Compartments are rectangular cuboids with the volume hh y h z where . . , N , k = 1, 2, . . . , K , be the number of molecules of the ith chemical species in the kth compartment at time t. Then, Z k (t) is an N -dimensional column vector with each component representing the number of molecules of the corresponding species in the kth compartment at time t. We define which is a K N-dimensional column vector and T denotes the transpose of a vector. We assume that the chemical system is subject to M chemical reactions with ζ j , j = 1, 2, . . . , M, being the corresponding N -dimensional stoichiometric vector. Let ζ k j , j = 1, 2 . . . , M, k = 1, 2, . . . , K , be a K N-dimensional stoichiometric vector which gives a net molecule change during each occurrence of the jth reaction in the kth compartment. Let ν k −,i (resp. ν k +,i ), i = 1, 2, . . . , N , k = 1, 2, . . . , K be a K N-dimensional stoichiometric vector which gives a net molecule change during diffusion of the ith species from the kth compartment to the (k − 1)th (resp. (k + 1)th) compartment. Let be the propensity function of the jth chemical reaction in the kth compartment, i.e. λ k j (Z k (t)) dt is the probability that the jth reaction occurs in the kth compartment during the time [t, t + dt) given that the current state at time t is Z k (t). We denote by R k j (t), j = 1, 2 . . . , M, k = 1, 2, . . . , K , a random process which counts the number of times the jth reaction occurs in the kth compartment up to time t. Then, where Y k j are independent unit Poisson processes. We define R k −,i (t) (resp. R k +,i (t)), i = 1, 2, . . . , N , k = 2, 3, ..., K (resp. k = 1, 2, ..., K − 1), random processes counting the numbers of times that one molecule of the ith species in the kth compartment diffuses to the (k − 1)th compartment (resp. to the (k + 1)th compartment) up to time t. Then, where Y k ±,i are independent unit Poisson processes. The governing equation for the state vector Z(t) is When the propensities are large (Kurtz 1978), the counting processes in Eqs.
(1)-(2) can be approximated as where W k j and W k ±,i are standard Brownian motions. Using ν k +,i = −ν k+1 −,i for k = 1, 2, . . . , K − 1 and changing the index (k + 1) → k in the last term of Eq. (3), governing Eq. (3) can be approximated by the following SDE (Gillespie 2000;Kurtz 1978) Since W k −,i and W k−1 +,i terms always appear together in Eq. (4), and since the sum of independent normal random variables is normally distributed, Eq. (4) can be rewritten as where W k−1 i is a standard Brownian motion. Let V h = hh y h z be the volume of each compartment, and define c(t) = Z(t)/V h as a concentration vector for species at time t. Define where c k i (t) = Z k i (t)/V h . The second part of Eq. (6) is consistent with the discretized Langevin scheme for a diffusion equation, as studied in Alexander et al. (2002). We rewrite Eq. (6) using the fact that reaction happens among species in the same compartment and that diffusion occurs between neighbouring compartments. Differentiating Eq. (6), the concentration of the chemical species in the kth compartment satisfies where W k (t) are N × N diagonal matrices with W k i (t) on its diagonal for i = 1, 2, . . . , N , k = 1, 2, . . . , K − 1 and χ {·} is an indicator function. In Eq. (7), ζ j is an N -dimensional stoichiometric vector of the jth reaction for j = 1, 2, . . . , M, and D is a N × N diagonal matrix which has diffusion constants of individual species on its diagonal, i.e.
x +Δx is normally distributed with zero mean and variance ΔxΔt. Matrices ξ (x, t) are diagonal N × N matrices where diagonal entries are independent spatio-temporal white noise processes. Then, Eq. (7) is a solution of a discretized version of a SPDE in space which can be formally written in the following form where c(x, t) is a spatio-temporal concentration related to c k (t) by The reaction termλ j : Note that Eqs. (6)-(7) are discretized versions of Eq. (8), but the compartment-based model in (3) breaks down as h → 0 as discussed in Section 2.2 of Engblom et al. (2009). The SPDE in Eq. (8) is consistent to the ones in the previous work (see Equation  1 in Kim et al. 2017 and Equation 3.24 in Dogan and Allen 2011). For more details, see derivations of the SPDE for diffusion in Section 3.1 of Dogan and Allen (2011)

Multiscale Algorithms Combining Compartment-Based Models with SPDEs
In this section, we present a multiscale approach which uses both SPDEs and Markov chain models. We develop two algorithms, denoted Scheme 1 and Scheme 2 in what follows, which are applied to illustrative examples in Sects. 4, 5 and 6. Considering the same set up as in Sect. 2, we study a system of N chemically reacting species The main goal of this paper is to replace the Markov chain description in a part of the computational domain by the corresponding SPDEs. Let us consider that we use the SPDE in Eq.
i.e. we consider that the first K s compartments are described by a suitable discretization of the SPDE in Eq. (8) (see Fig. 1b). We only use the Markov chain model for the remaining K m = K − K s compartments, i.e. in subdomain In this section, we develop an appropriate boundary condition on the interface I between Ω s and Ω m .
In order to design the numerical scheme, the SPDE in Eq. (8) needs to be appropriately discretized. We denote by Δx the mesh size used in the discretization of the SPDE. There are two important cases: (i) Δx > h and (ii) Δx ≤ h. In this section, we focus on case (ii), because we are interested in coupling the SPDE in Eq. (8) with Markov chain models. The case (i) is important when one uses discretized SPDEs to design efficient multiscale schemes, but this introduces additional discretization errors. We will discuss case (i) in Sect. 7. In Ω s , each compartment of size h is discretized into α grid points (α ∈ N) with each grid size equal to Δx. In the remaining part of the computational domain Ω m , the compartment-based model is used. The state of the system of the multiscale model is described by vectors The vector X K s α+k (t) for k = 1, 2, . . . , K m represents species numbers in C K s +k = We consider two different schemes to describe transfer of molecules near the interface I coupling discretized SPDEs and the Markov chain model, as shown in Fig. 2. Without loss of generality, both schemes are introduced for diffusion, because the description of reactions does not influence the transfer of molecules across the interface I . In Scheme 1, we assume that there is a virtual compartment, in Ω s , where the molecules are partially treated using a compartment-based approach. Such overlap (handshaking) regions are common in many multiscale methodologies, including coupling molecular dynamics with Brownian dynamics simulations (Erban 2014(Erban , 2016, Brownian dynamics with PDEs (Franz et al. 2013), or in atomistic to continuum coupling methods (Miller and Tadmor 2009). We define a state vector which is a (K s α + K m )N -dimensional column vector. Scheme 1 is described using the following evolution equation for state vector X(t): where the first sum on the right-hand side represents diffusion in Ω s (compare with Eq. (5) replacing h by Δx). Note that the indicator function χ is used to make sure the term inside the square root not being negative. Here the symbols ν k ±,i describe (K s α+K m )N -dimensional stoichiometric vectors. The second and third sums represent diffusion in the compartment-based region, Ω m , where Y k ±,i are independent unit Poisson processes (compare with Eq. (2)). The last two sums represent transition of a molecule from Ω m to Ω s and from Ω s to Ω m , respectively. A molecule in Ω m in the boundary compartment, C K s +1 , jumps to the SPDE domain with a rate D i /h 2 . A molecule which jumps is placed to one of the mesh points in the overlap compartment, C K s . To describe this process in Eq. (9), we have defined indicator functions where U ±,i (t) are independent random variables uniformly distributed on interval [I − h, I ] for each t and i. Stoichiometric vectors, η ±,i for = 1, 2, . . . , α, i = 1, 2, . . . , N , give changes due to the diffusion of the ith species between the th SPDE discretization point in C K s and the compartment C K s +1 across the interface I . Transition of a molecule from Ω s to Ω m is described by the last term of Eq. (9) using time-changed Poisson processes. A molecule, anywhere in the overlap compartment C K s , can be transferred with a rate D i /h 2 . The corresponding molecule is then randomly subtracted from one of α discretization grid points which are in C K s . Note that the molecular copy number, α j=1 X in the last term of Eq. (9) can be non-integer value due to the non-integer concentration in C K s . To prevent X (K s −1)α+ j i (s) being negative due to the molecular transfer from Ω s to Ω m , another indicator function is used in the last term of Eq. (9) to set the propensity as zero if the total molecular copy number in C K s is less than 1.
Scheme 2 is described in terms of two unknown parameters, denoted Ψ 1 and Ψ 2 , by the following evolution equation for the state vector X(t): The first three sums on the right hand side in Eq. (10) are identical to those in Eq. (9). The fourth and fifth sums describe molecular transfer between the last grid point in Ω s and the boundary compartment C K s +1 . A molecule in Ω m in the boundary compartment, C K s +1 , jumps to the last grid point of the SPDE domain with rate Ψ 1 D i /h 2 , and the transfer rate in the opposite direction is Ψ 2 D i /Δx 2 . Note that X K s α i (s) in the fifth term of Eq. (10) can be non-integer value due to the non-integer concentration in Ω s . To prevent X K s α i (s) being negative due to the molecular transfer from Ω s to Ω m , we use an indicator function to set the propensity as zero if the molecular copy number in the last grid point in Ω s is less than 1.
To determine parameters Ψ 1 and Ψ 2 of Scheme 2, we use the discretization of the PDE for diffusion using a finite volume approximation (Bernstein 2005). It gives the jump coefficient of the ith species from the jth compartment to the neighbouring j th compartment as D i /(h j |a j − a j |), where h j is the length of the jth compartment and a j and a j are the centres of the jth and j th compartments, respectively. Considering the size of the domain allowed for molecule transfer across the interface in Scheme 2, we set |a j − a j | = (Δx + h)/2. We take h j = Δx for the jump coefficient from Ω s to Ω m and h j = h for the jump coefficient from Ω m to Ω s . Then, we match the jump coefficients to the rate constants for jump across the interface given in Eq. (10) to derive the following formula for the parameters of Scheme 2 The multiscale algorithm for Scheme 1 for the case of diffusion only is given in Table 1.
We denote a propensity of diffusion of the ith species in the This definition also includes the propensity of a diffusive jump of the ith species from the Markov chain domain, given as a 1 −,i (t). We denote a propensity of diffusive jump of the ith species from the SPDE domain by Table 1 Pseudocode for the multiscale reaction-diffusion algorithm with Scheme 1 applied to simulation of diffusion Then, we define total propensity in Ω m Total propensity a 0 is used in steps [A] and [B] in the pseudocode in Table 1 to select time when the next event occurs in Ω m . The pseudocode denotes the time of the next update in each subdomain as t s and t m , and the current time as t. In step [B], we update the compartment-based part of the system. In step [C], we update the SPDE part of the system by where ζ k−1 i are independent normally distributed random numbers with zero mean and unit variance.

Application: Static Boundary
In this section, we apply the multiscale approach to examples in which we know a priori the position of the boundary I between Ω s and Ω m . Generalization to a more complicated case with a moving boundary is presented in Sect. 5.

A Morphogen Gradient Model
We consider a morphogen gradient model in It consists of one chemical species S, i.e. Z k (t) is a scalar describing the number of molecules of S in C k . The state of the Markov chain model is described by the K -dimensional column vector Z(t) = Z 1 (t), Z 2 (t), . . . , Z K (t) T . Morphogen is subject to diffusion which is described by Eq.
(2). There are also two reactions in the system. Morphogen, S, is produced in the first compartment with rate J , i.e. the propensity is λ 1 1 Z 1 = J . Morphogen degrades everywhere with rate δ, i.e. with propensity λ k 2 Z k = δ Z k for k = 1, 2, . . . , K . In all stochastic simulations of the morphogen gradient model, we assume that 500 morphogen molecules are initially uniformly distributed in the half of the domain, Ω s = [0, L/2] × [0, h y ] × [0, h z ]. The parameters are given in Table 2.
where D is the diffusion constant of morphogen S. We apply the multiscale approach using both schemes developed in Sect. 3. Since the morphogen is produced at the left end, the morphogen has a decreasing gradient as it goes towards L. Therefore, we split the spatial domain in half and set the left half as Ω s and the right half as Ω m , i.e. I = L/2. The (K s α+ K m )-dimensional state vector of the multiscale model is denoted X(t) = X 1 (t), X 2 (t), . . . , X K s α+K m (t) T . Note that the morphogen molecules are only produced in the first discretization mesh point with size Δx in Ω s . In Fig. 3, we simulate the morphogen gradient model using Scheme 1 of the multiscale algorithm. We calculate 10 4 realizations of the sample paths of the stochastic process and present mean and standard deviations of the morphogen numbers in Ω at different times, t = 0, 2, 5, 20 s. Morphogen numbers in α grid points of Ω s are summed so that they can be compared to the numbers in the underlying Markov chain model. We compare the results with mean and standard deviations of the morphogen numbers which we calculate analytically using matrix analysis for reaction-diffusion Markov chain models (Gadgil et al. 2005;Kang et al. 2012b). In Fig. 3, morphogen numbers in Ω s (resp. in Ω m ) are expressed as green bars (resp. blue bars). Error bars represent one standard deviations from the mean number of morphogens in each compartment. Mean and standard deviations of the morphogen numbers from the analytic solution are drawn as a red line and blue dotted lines. The results using the multiscale algorithm match perfectly to the ones from the exact solution.
In Fig. 4, we present relative errors of the means and standard deviations of the number of molecules between the Markov chain model and multiscale model. The analytic solution is used for the statistics of the Markov chain model, and both schemes where E[·] and σ [·] represent a mean and standard deviation. In Fig. 4a, red and green lines represent e m (k) and e v (k) at time t = 50 s using Scheme 1, respectively, and blue and purple lines are for Scheme 2. We observe that the relative errors in Eq. (13) are less than 4 % in the entire simulation domain. In Fig. 4b, we compare the maximum absolute values of the relative errors defined in Eq. (13) with α = 1, 5, 10, 25 and fixed compartment size h where α = h/Δx. In both schemes, the relative errors are in a range of less than 4 % except for the case when α = 25 with Scheme 2. The relative errors in the mean and standard deviation become significantly larger when we apply the multiscale algorithm using Scheme 2 with α = 25. In this case, the mean in C K s +1 gets larger than the mean in C K s which shows a bias in the method for larger values of α (the exact mean number of molecules decreases along the x-axis). We provide an explanation of this phenomenon in the next section.

A Diffusion Model with Two Compartments
In Sect. 4.1, we have observed that the error of Scheme 2 increases when we decrease the ratio of the numerical discretization in Ω s and the compartment size. When α = 25, the mean number of molecules of the morphogen does not have a decreasing gradient across the interface I in Scheme 2. To investigate this numerical error, we consider diffusion in a theoretical two-compartment model. It is similar to the one in Sect. 4.1, but we set J = δ = 0 so that there is no flux or degradation of the morphogen. We use L = 2h and Then, each region consists of one compartment, K s = K m = 1, and X(t) is an (α + 1)-dimensional vector. In Fig. 5a, we present simulation results of the two-compartment model using Scheme 1 (red line) and Scheme 2 (green line) with α = 10, 20, 30, 40, 50 and compare them to the simulation result of the Markov chain model using the Gillespie SSA (purple line). The Markov chain model has α + 1 numerical grid points where the first α ones are with size Δx = h/α and the last one with size h. Diffusion of molecules is simulated by jumps from grid points to their nearest neighbours, i.e. the numerical meshes in the Markov chain model are coupled by diffusion in the same way as it is done in Scheme 2. Applying both multiscale algorithms and the Gillespie SSA, we compare the mean morphogen numbers in the second compartment computed from 100 realizations of simulation. Using 50 molecules in total, the exact value of the mean numbers of molecules in C k , k = 1, 2, is 25. Notice that Scheme 1 and the Gillespie SSA with two mesh sizes correctly approximate the means. However, Scheme 2 overestimates the mean morphogen number in C 2 as α gets large. To understand where the numerical error arises, we also simulate Scheme 2 without the noise term in the SPDEs (marked as a blue line in Fig. 5a), i.e. we remove the term with a square root in Eq. (10).
In Fig. 5a, we observe that the mean morphogen number in C 2 , E[X α+1 ], is underestimated as α increases when we use Scheme 2 without noise term in the SPDEs. Note that X , = 1, 2, . . . , α, always have non-negative integer values due to no noise term in modified Eq. (10). The molecular transfer from Ω s to Ω m occurs when X α ≥ 1. However, the frequency of this transfer is not sufficient as α gets large, which lowers X α+1 . On the other hand, with noise terms included in Eqs. (9) and (10), there are more chances that X < 0 for some = 1, 2, . . . , α due to large fluctuations with a small number of molecules as α gets large. Then, it is more frequent that X α ≥ 1 due to the fact that α+1 =1 X = 50 and X < 0 for some = 1, 2, . . . , α. More frequent molecular transfer from Ω s to Ω m causes overestimation of the mean morphogen number in C 2 in Scheme 2.
In Fig. 5b, we compare distributions of the morphogen numbers when α = 10 and 50. The distributions are computed from 1000 realizations of simulation when t = 0.01 s and 50 s. Each distribution is computed over all X , = 1, 2, . . . , α so that we can display an overall range of the morphogen number in each discretization of Ω s . Each X is normalized by Δx so that the distributions can be compared for different α's. The normalized mean morphogen number (density) in Ω s decreases significantly in both cases with α = 10 and 50 as time evolves. On the other hand, the variance of the morphogen density is much greater for α = 50 than for α = 10 at t = 0.01 s due to the lower morphogen number in each discretization of Ω s . Therefore, we conclude that the error in Scheme 2 strongly depends on the size of fluctuations close to the interface. On the other hand, the molecular transfer from Ω s to Ω m is decided by α =1 X in Scheme 1. This makes Scheme 1 more robust than Scheme 2 for large values of α.

Application: Moving Boundary
In some applications (Robinson et al. 2014), it is difficult to decide a position of the interface I a priori. In this section, we extend the presented algorithm to the case when the location of the interface I (t) between Ω s and Ω m moves in time, based on the number of molecules in each part of the domain. The multiscale approach with the adaptive interface is applied to the example introduced in Sect. 4. The adaptive algorithm is described in Table 3. Following (Robinson et al. 2014), we introduce two thresholds denoted Q upper and Q lower (Q upper ≥ Q lower ), and one integer parameter n c . We initialize the position of the interface I (0) = 0 in step [A'], i.e. we initially model the whole domain using the detailed compartment-based approach. We run the original Scheme 1 until time n c Δt. We check whether the interface I (t) should be moved in step [C']. If the number of molecules in the compartment next to the interface in Ω s is smaller than Q lower , a compartment-based model is used in that region. On the other hand, if the number of molecules in the boundary compartment next to interface I (t) in Ω m is larger than threshold Q upper , the corresponding compartment is transferred to the SPDE region where the molecules are redistributed uniformly in α grid points. Due to the uniform redistribution of the molecules, rapid changing of the interface I (t) introduces more errors. Note that in Scheme 1 with a fixed boundary, one molecule has been chosen randomly from α discretizations of C K s in Ω s and transferred to C K s +1 in Ω m . Similarly, we have taken one molecule from C K s +1 and transferred the molecule to the randomly chosen SPDE mesh point in C K s . However, in the adaptive algorithm, we modify the setting of Scheme 1 so that a molecule is taken uniformly from the entire region of C K s and transferred to C K s +1 , i.e. 1/α molecule is subtracted in all α SPDE grid points of C K s . Similarly when the molecule is transferred from C K s +1 to C K s , 1/α molecule is added in all α grid points of C K s . Without this modification of the setting in Scheme 1, the appropriate level of the morphogen gradient is not formed in Fig. 6.

The adaptive algorithm [A']-[D'
] is applied to the morphogen gradient model introduced in Sect. 4, and the results are presented in Fig. 6. We use Q lower = 15, Q upper = 25 and n c = 10. Other parameters are given in Table 2. Our initial condition is X k (0) = 0, for k = 1, 2, . . . , K s α + K m , i.e. the system starts with no molecules and the gradient is formed during the simulation. In Fig. 6, one realization of the algorithm in Table 3 at different times t = 0.5, 2, 10, 40 s is presented. The green and blue bars represent the numbers of molecules in the corresponding compartments in Ω s and Ω m , respectively. The blue dotted line represents interface I (t) between two regions, and the red circles are the mean numbers of molecules obtained from the analytic solution of the stochastic model. Our results show that the boundary between two regions is moving to the right in time as the molecule numbers increase due to the production on the left.
In Fig. 7a and b, we simulate the adaptive algorithm with fixed thresholds for a range of values of n c = 1, 10, 10 2 , 10 3 , 10 4 , which are the numbers of time steps to Table 3 Pseudocode for the adaptive multiscale reaction-diffusion algorithm with Scheme 1 applied to simulation of diffusion  Table 1 until time k(n c Δt)

[D'] Repeat steps [B']-[C']
until the simulation ends Fig. 6 Comparison between one realization of the number of morphogens using Scheme 1 with a moving interface, given in Table 3 (green bars and Fig. 7b and c, we observe that the maximum absolute values of the relative errors increase as the number of time steps, n c , or the size of the threshold window, Q upper − Q lower , gets smaller. This is because the small size of the number of time steps (or the threshold window) makes the interface location change frequently, which causes additional errors. On the other hand, Fig. 7a and d do not show similar pattern since large size of the threshold window (Q upper − Q lower = 10) and the number of time steps (n c = 10 3 ) prevent frequent movement of the interface location. Overall, Scheme 2 has slightly smaller errors than Scheme 1. In Fig. 7, the maximum absolute values of the relative errors are calculated using 10 4 realizations of simulation using Scheme 1 or 2 for each value of n c and for each set of values of (Q lower , Q upper ) and using the analytic solution of the Markov chain model.

Applications: Multiple Species
In this section, we illustrate the applicability of the presented multiscale approaches to chemical systems with multiple species. Since different chemical species can have very different molecular distributions in the computational domain, the partition of the computational domain into subdomains Ω s and Ω m can be species dependent. We use the pom1p gradient model from Saunders et al. (2012) to illustrate a multiscale approach, where each species has a different partition into Ω s and Ω m depending on its molecular distribution. The model consists of two species, slow-diffusing pom1p clusters, denoted S 1 , and fast-diffusing pom1p particles, denoted S 2 . We use pseudo 1-dimensional domain Ω as in Fig. 1, where L = 14 μm, which is divided into K = 40 compartments, C k , k = 1, 2, . . . , K . Both S 1 and S 2 are produced in the whole computational domain with space-dependent rates given by Saunders et al. (2012), i.e. with propensities λ k j (Z k ) = a j exp −a 6 k − K + 1 2 where j = 1, 2, k = 1, 2, . . . , K , and a 1 , a 2 and a 6 are constants given in Table 4. In addition to production, species S 1 and S 2 are subject to the following reactions which take place in the whole domain with the corresponding propensities given by where k = 1, 2, . . . , K , and a 3 , a 4 and a 5 are constants given in Table 4. In Fig. 8, we present an illustrative simulation of pom1p gradient model. We plot spatial distributions of S 1 and S 2 at times t = 50 s and t = 1000 s. We observe that the spatial distribution of S 1 contains a region with high abundance of molecules in the centre of the computational domain. The chemical species S 2 has low copy numbers in the entire domain. Therefore, we introduce the SPDE region in the middle of the domain by (note that we fix K = 40 in this example) where the coarse graining is only applicable to S 1 in Ω s . In particular, we have introduced two interfaces, I 1 and I 2 between Ω s and Ω m . Diffusion of chemical species S 1 is simulated using the algorithm in Table 1. Similarly, production of S 1 is implemented using the SPDE and Markov chain model in Ω s and Ω m , respectively, as it is done in Eq. (9). The chemical species S 2 is simulated by the Markov chain model in the entire domain, because the average number of molecules of S 2 is relatively low. In particular, diffusion, production and degradation of S 2 are implemented as in the underlying Markov chain model. The only complications are reactions because they include both species S 1 and S 2 , which are in Ω s described by different modelling approaches. We treat these reactions as time-changed Poisson processes in both subdomains Ω m and Ω s . Discretizing each compartment, C k , k = 11, 12, . . . , 30, into α grid points, the state of S 1 variable is described by vector, X(t) = (X 1 , X 2 , . . . , X 20(α+1) ) where X 1 , X 2 , . . . , X 10 (resp. X 20α+11 , X 20α+12 , . . . , X 20(α+1) ) are the numbers of molecules of S 1 in the left (resp. right) part of Ω m . The values of SPDE description in compartment C k , k = 11, 12, . . . , 30, are given by X 10+(k−11)α+ , = 1, 2, . . . , α. The state of S 2 variable is described by vector, where Y k is the number of molecules of S 2 in compartment C k , k = 1, 2, . . . , K . The propensity of the first reaction in (14) of the multiscale model is given by The propensity of the second reaction in (14) of the multiscale model is given by We simulate reactions in (14) as time-changed Poisson processes with propensities in Eqs. (15)-(16). If the first of these reactions occurs in C k , k = 11, 12, . . . , 30, we subtract 1/α from each X 10+(k−11)α+ , = 1, 2, . . . , α, and we add one to Y k . If the second reaction in (14) occurs in C k , k = 11, 12, . . . , 30, we add 1/α to each X 10+(k−11)α+ , = 1, 2, . . . , α, and we subtract one from Y k . Note that the conversion of S 1 in C k , k = 11, 12, . . . , 30, is applied equally to the entire α grid points of C k rather than to one randomly chosen grid point in C k as it is done (consistently with Eq. (9)) for diffusion across the interfaces. In Fig. 8, green bars and blue bars represent the mean numbers of molecules of the pom1p clusters and particles in Ω s and Ω m using the multiscale algorithm with Scheme 1. Error bars represent one standard deviation from the mean in the multiscale approach. Red lines and blue dotted lines are the mean numbers and their standard deviations from the means computed by the Gillespie SSA simulating the compartment-based approach in the entire domain. Both statistics using the compartment-based approach and the multiscale algorithm are computed by averaging over 10 4 realizations of the simulations for each case.

Discussion
A Markov chain model (compartment-based model) has been widely used to describe the discrete nature of the molecular copy numbers and inherent stochasticity in reaction-diffusion systems (Erban and Chapman 2019), but it can be computationally intensive. A possible approach to increase efficiency of simulations is to approximate a part of the model by some coarse-grained methods. In this paper, we have introduced two multiscale algorithms coupling the SPDEs and the Markov chain model, which provide good approximations to the solutions obtained by the Markov chain model applied in the entire spatial domain. Two coupling methods of the Markov chain model and the SPDEs across the interface have been studied. In this section, we compare the presented approach with methods in the literature.
Several Langevin formulations have been introduced to model fluctuating hydrodynamics for chemically reactive species (Bhattacharjee et al. 2015) and stochastic reaction-diffusion systems (Kalantzis 2009;Ghosh et al. 2015). In particular, the spatial chemical Langevin equation was applied to the Gray-Scott model, and its pattern formation was compared to the ones obtained by the reaction-diffusion master equation and PDEs (Ghosh et al. 2015). The spatial chemical Langevin equation consists of a system of stochastic differential equations, and it corresponds to Eq. (4) in Sect. 2. On the other hand, several approaches using SPDEs (Atzberger 2010;Dogan and Allen 2011;Alexander et al. 2002Alexander et al. , 2005 have been introduced to model stochastic reaction-diffusion systems. In Atzberger (2010), the SPDE was derived for reaction-diffusion systems, and discretization of PDEs and stochastic fields was discussed. Unlike Eq. (5), the stochastic fields in the discretized SPDEs account for fluctuations due to diffusion but not for reaction. In , the SPDE for reaction-diffusion systems was derived which is consistent with Eq. (8). In their formulation, diffusion was implemented by the SPDE while the reaction was simulated using the exact or modified SSA.
In Yates and Flegg (2015), two hybrid algorithms are suggested for coupling a compartment-based model and a PDE model when the size of the PDE discretization is less than or equal to the compartment size. Both algorithms extend the PDE approach to the systems with low copy numbers of molecules in a part of the computational domain. The first algorithm considers the PDE solution as the probability density to find a particle within the region and is applied to both cases of low and high copy numbers of molecules in the PDE region. The second algorithm is a simplified and more efficient version of the first one when the PDE region involves the high copy number of molecules. Like in this paper, both algorithms implement a pseudocompartment with size h in the PDE region where h represents the compartment size. The second algorithm in Yates and Flegg (2015) is similar to Scheme 1 if a discretized version of SPDEs replaces the PDEs. However, the interface between the two modelling regimes is assumed to be fixed in Yates and Flegg (2015). In Harrison and Yates (2016), a hybrid algorithm is introduced coupling a compartment-based model and PDEs where the size of the PDE discretization is much finer than the compartment size. In the model, an overlap region is defined with two interfaces (corresponding to the pseudo-compartment in Scheme 1) where both modelling regimes are valid, and both cases with fixed and adaptive interfaces are considered. Unlike our pseudocompartment in Scheme 1, the overlap region can contain multiple compartments if needed. On one interface between the compartment-based model and the overlap region, the population of the PDE solution on the interface is matched to the average of the population in the neighbouring compartments. On the other interface between the PDE region and the overlap region, flux on the interface was matched. The hybrid algorithm in Harrison and Yates (2016) approximates the mean population numbers in the compartment-based model if it was possible to apply it over the entire spatial region. The use of the overlap region allows matching the variance between two models in the compartment-based region when the fixed interface is used. On the other hand, the goal of Scheme 1 and Scheme 2 is to approximate the compartment-based model by employing the discretized version of SPDEs in the region with high molecules. Therefore, we can match both the mean and variance of the population numbers computed by our multiscale algorithms to the results obtained by the compartment-based model in the whole spatial domain. This is done for both cases with a fixed or adaptive boundary. Unlike the previous approaches in Yates and Flegg (2015); Harrison and Yates (2016), the presented multiscale algorithms can apply to systems with multiple species as it is shown in Sect. 6 where each species has a different partition of the spatial domain into subdomains where different models are used, depending on the spatial distribution of molecules of each species. In Spill et al. (2015), a hybrid algorithm is presented using a compartment-based model and PDEs, where the size of the compartment and numerical discretization for the PDE model is equal.
In this paper, we have discussed the case when the mesh size of the numerical discretization of the SPDEs is smaller (or equal) than the compartment size in the Markov chain model (h ≥ Δx). This case is useful when we add inherent stochasticity in the PDE model where a fine spatial resolution of the PDE solution is required to describe the solution of the SPDE. This case was also discussed in other hybrid algorithms coupling the compartment-based model and the macroscopic PDEs (Yates and Flegg 2015;Harrison and Yates 2016). The other case, h < Δx, discussed, for example, in the hybrid algorithm coupling a random walk on a lattice and the PDE model (Flekkoy et al. 2001), is helpful when the PDE or SPDE model is used as a coarse-grained approximation of the compartment-based model. Such approximation can be used in the region where spatial concentration gradients are not large, so they do not require a fine resolution in space. Although we have focused on the case h ≥ Δx, the presented approach can be extended to h < Δx as well. In fact, if h = Δx, both Scheme 1 and Scheme 2 will be the same. If h < Δx, we may be able to consider an overlap region (like a pseudo-compartment) in the compartment-based region to extend Scheme 1. The presented SPDE-based approach provides a bridge between the stochastic approach (using the Markov chain compartment-based model) and the deterministic approach (using the macroscopic PDEs) by incorporating a discretized version of SPDEs. The SPDEs can be utilized to build other hybrid models, for example, by coupling them with macroscopic PDEs. Then, some approaches used in the hybrid algorithms coupling the compartment-based model with the PDEs (Yates and Flegg 2015;Harrison and Yates 2016;Spill et al. 2015;Smith and Yates 2018) will naturally apply to the case with the SPDEs.