A hybrid method for micro-mesoscopic stochastic simulation of reaction-diffusion systems

Mathematical Biosciences 312 (2019) 23–32

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

A hybrid method for micro-mesoscopic stochastic simulation of reactiondiffusion systems

T

Alireza Sayyidmousavi , Katrin Rohlf, Silvana Ilie ⁎

Department of Mathematics, Ryerson University, 350 Victoria St, M5B 2K3 Toronto, Canada

ARTICLE INFO

ABSTRACT

Keywords: Hybrid method Stochastic simulation Multi-Particle Collision dynamics Inhomogeneous Stochastic Simulation Algorithm Reaction-diffusion systems

The present paper introduces a new micro–meso hybrid algorithm based on the Ghost Cell Method concept in which the microscopic subdomain is governed by the Reactive Multi-Particle Collision (RMPC) dynamics. The mesoscopic subdomain is modeled using the Reaction-Diffusion Master Equation (RDME). The RDME is solved by means of the Inhomogeneous Stochastic Simulation Algorithm. No hybrid algorithm has hitherto used the RMPC dynamics for modeling reactions and the trajectories of each individual particle. The RMPC is faster than other molecular based methods and has the advantage of conserving mass, energy and momentum in the collision and free streaming steps. The new algorithm is tested on three reaction-diffusion systems. In all the systems studied, very good agreement with the deterministic solutions of the corresponding differential equations is obtained. In addition, it has been shown that proper discretization of the computational domain results in significant speed-ups in comparison with the full RMPC algorithm.

1. Introduction Stochastic modeling and simulation have been successfully used to study many biological and physical systems. This is particularly important when there are small numbers of particles of a particular chemical species in a biochemical system leading to inherent molecular fluctuations, in which case deterministic models are often inaccurate if even applicable. A number of stochastic modeling algorithms have so far been proposed and used by researchers. Stochastic Simulation Algorithm (SSA) is one of the most frequently utilized examples of such stochastic schemes which was originally proposed by Gillespie [1,2]. Later, other researchers [3–5] advanced the Stochastic Simulation Algorithm to spatial systems leading to the Inhomogeneous Stochastic Simulation Algorithm (ISSA), for solving the Reaction-Diffusion Master Equation. A recent review on RDME has been conducted by Smith and Grima [6]. Solution of the RDME using ISSA involves the discretization of the problem into a number of subvolumes each being treated as a homogeneous system. Diffusion is modeled as discrete jumps between neighboring subvolumes. In reality, particles can freely move in a system and collide with one another resulting in the change of their velocities and identities [7]. Elf and Ehrenberg [8] introduced the Next Subvolume Method as one of the most common methods to implement the ISSA. Both the accuracy and speed of the ISSA have been improved

⁎

by the development of hybrid and adaptive algorithms [9–11]. The ISSA has also been enhanced by some researchers to include reactions with delays [12,13]. Nevertheless, the fact that ISSA is incapable of keeping track of single particles, renders it applicable mainly when there is no need for tracking particles and detailed information on the interaction between the particles. A more detailed approach for modeling microscopic systems are off-lattice Molecular Dynamics models [14] where individual particle motion can be traced. This method, however, is computationally expensive when trying to simulate a chemical reaction in a solvent [15]. Brownian Dynamics (BD) provides a more efficient alternative for simulating diffusion as used in Smoldyn [16], as do Green's Function Reaction Dynamics algorithms [17,18]. A more recent less detailed particle-based method, named MultiParticle Collision (MPC) dynamics, is an off-lattice algorithm introduced by Malevanets and Kapral [19,20]. The MPC dynamics was later extended to include reactions giving rise to the Reactive MultiParticle Collision (RMPC) dynamics [21,22]. In many biological systems, accuracy and microscopic details are required only in parts of the computational domain while the rest of the domain can be modeled on a more coarsely grained level. Modeling such systems calls for the development of hybrid algorithms, which can simulate different parts of the computational domain at different levels of detail. The most important benefit of such hybrid approaches is to speed up simulations and thereby save considerable computational

Corresponding author. E-mail address: [email protected] (A. Sayyidmousavi).

https://doi.org/10.1016/j.mbs.2019.04.001 Received 24 August 2018; Received in revised form 13 April 2019; Accepted 14 April 2019 Available online 15 April 2019 0025-5564/ © 2019 Elsevier Inc. All rights reserved.

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

time. While the term hybrid encompasses methods that couple simulations at different scales, the focus of the present paper is on coupling microscopic dynamics to mesoscopic dynamics. Micro–meso hybrid models are required when fluctuations do not need to be accounted for through the whole domain, but where specific particle details are not necessary only in some parts of the domain [23]. The importance of considering fluctuations necessitates the use of stochastic algorithms at both micro and meso levels. In general, hybrid simulations can be performed using the adjacent-space technique, or the overlapping space. In the latter method, both algorithms are used in the 'overlapped' region of space, but applied to different biochemical species in the system. Several micro–meso hybrid methods have so far been developed. Flegg et al. [24] proposed the so-called Two Regime Method. The domain is partitioned into meso and micro regimes separated by an interface. In the micro regime the particle paths are modeled using Brownian dynamics, while the meso regime is simulated using an eventdriven compartment-based approach. A probability theory is used to model the flux over the interface so that no error occurs if the net flux is zero. The Two Regime Method was then extended by Robinson et al. [25] to model an adaptive interface allowing for a dynamic selection of subdomains where the more refined Brownian dynamics simulation is needed. The Two Regime Method has also been extended to higher dimensions by Flegg et al. [26] where the authors modeled a square domain with interfaces which are either horizontal or vertical. Another micro–meso hybrid model, called the Compartment-Placement-Method, was presented by Hellander et al. [27]. In this model, the problem domain is discretized into subcells. First, the particles in the microscopic portion of the domain are frozen and the particles in the mesoscopic portion are simulated using the ISSA. Next, the microscopic particles advance in accordance with Green's Function Reaction Dynamics as an alternative to Brownian dynamics. The Two Regime Method and Compartment-Placement-Method as discussed by Flegg et al. [28] do not converge in the limit of small time steps and also when compartment size is kept fixed. As a result, they proposed another technique named the Ghost Cell Method (GCM). Similar to other micro–meso hybrid models, in the GCM, the domain is divided into micro–meso subdomains simulated using the Brownian dynamics and a compartment based method, respectively. The two subdomains are separated by a reflective interface and a ghost cell is constructed in the microscopic domain to allow the movement of particles in between the two regimes. The present paper proposes a new micro–meso adjacent-space hybrid method based on the GCM concept with the microscopic subdomain modeled using the RMPC dynamics. The mesoscopic subdomain is modeled using the RDME. No hybrid algorithm of this type has hitherto employed the RMPC dynamics for modeling reactions and the trajectories of each individual particle. The MPC uses two steps in its algorithm: free-streaming, and collisions. Unlike for lattice-based methods, positions and velocities of the particles are continuous in space and time, and not constrained to lie on an underlying lattice structure. To simulate local collisions, the domain is discretized into smaller subvolumes. The collision is a result of a random rotation operator acting on the velocities of the particles in the subvolume in such a way that mass, momentum and energy are conserved in the subvolume. With respect to kinetic theory literature, MPC would be classified as being a mesoscale model, while it would be called microscopic in so far as stochastic chemical reaction literature is concerned. The MPC collision step is less detailed than that of Lattice-Boltzmann [29] or Dissipative Paticle Dynamics [30] methods, and thus also more efficient. Ihle and Kroll [31–33] and Noguchi and Gompper [34] are among the researchers who further analyzed the MPC from a theoretical viewpoint. The MPC dynamics has been utilized to simulate different systems with hydrodynamic flows [35], polymer systems [36,37], colloid [38] and self-propelled particles [39,40]. Using the MPC, Tucci and Kapral [21,22] studied diffusion-influenced reaction kinetics leading to

the development of the so called Reactive MPC (RMPC). Rohlf et al. [7,41] advanced the RMPC dynamics further to simulate reactive mechanisms in which the total number of particles may change throughout the simulation. In a new study, Strehl and Rohlf [15] augmented the RMPC dynamics through the development of theoretical expressions for reactant-dependent diffusion control by taking into account speciesdependent diffusion coefficients. When compared to Molecular Dynamics simulations based on Brownian Dynamics, the RMPC is observed to have two major benefits. First, the number of random numbers generated in BD is of the order of N. In RMPC, however, it is of the order of NV × Nsp, where N, NV and Nsp denote the total number of particles, the number of subvolumes and the number of species types in the system, respectively. In other words, fewer random numbers are generated in the RMPC compared to BD. Second, the conservation of mass, momentum and energy is maintained in the collision and free streaming steps in RMPC which results in correct Navier–Stokes hydrodynamics. However, in BD, conservation of properties cannot be accounted for without the inclusion of hydrodynamic interactions. When it comes to modeling chemical reactions, the simulation of zeroth and first order reactions in BD involves the generation of Poisson and exponentially distributed random numbers. Modeling second order reactions in BD is based on the assumption that molecules can react only if they are within a specific distance. This assumption requires that the particles be initially positioned in the system such that the products, for instance in the case of a reversible reaction, do not unphysically react [42]. In RMPC, however, uniform random numbers are generated to decide on local reactions, which is less computationally costly compared to generating Poisson and exponentially distributed random numbers. Furthermore, the simulated diffusion in an RMPC simulation is generally calculated a posteriori while the current paper uses a strategy recently proposed by the authors in [43] that maintains a constant diffusion. Most reaction-diffusion models assume this to be true, and often treat the constant diffusion coefficient as an input to the model rather than measured a posteriori. This paper is organized as follows: The RMPC dynamics is described in Section 2. Section 3 describes the details of the compartment-based modeling using the ISSA. The algorithm of RMPC–ISSA hybrid model is presented in Section 4. In Section 5, the proposed hybrid approach is tested on three reactiondiffusion systems. Finally, the discussions are given in Section 6. 2. Reactive Multi-Particle Collision (RMPC) The RMPC comprises three main steps occurring at discrete time steps: particle collisions, reactions, and free-streaming. The system contains N particles in a finite volume V that is subdivided into NV subvolumes often also referred to as cells. Initial positions are assigned to the particles so as to correspond to a given initial distribution. Each of the three velocity components of the particles is sampled from a Maxwell–Boltzmann velocity distribution that has a zero mean and variance KBT/m. Here, KB T is the Boltzmann constant, T the temperature, and m the particle mass. A rotation operator ^ randomly chosen from a set of angles is assigned to each subvolume [44]:

lx 2 + (1 lx 2) c ^ = lx l y (1 c ) + lz s l x l z (1

c )

l x l y (1

lz s ly

l y lz (1

c ) + lx s

l y + (1

ly s

c )

2

2) c

lx lz (1

c ) + ly s

l y lz (1

c )

l z 2 + (1

lx s lz 2) c (1)

1 l y = sin 1 where l x = cos and lz = , Φ and θ represent random numbers on the interval [0,2π] and [−1,1], respectively and ψ is the collision angle controlling the diffusion. For brevity, the short form s = sin and c = cos is used. The velocity of the mass center of a subvolume is given as 2,

24

2

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

V =

1 N

diffusion systems. Reactions lead to changes in density of chemical species that thus affect the value of the diffusion coefficient locally, leading to non-constant, or rather, spatially dependent diffusion coefficients. If the diffusion coefficient is meant to be constant throughout the system volume for the duration of the simulation, one may have to compute the collision angle in a local subvolume at each time step so that Eq. (9) remains valid throughout the simulation. The second issue is that Eq. (9) is valid only for systems in which the particle density is uniform in the simulation volume, that is, when particles are distributed uniformly in the system. This latter issue can be resolved through introduction of a uniformly distributed set of bath particles that then act to provide a large enough local density so that Eq. (9) holds for all subvolumes regardless of any spatial variations of density as a result of inhomogeneous initial conditions or reaction depletion events. The bath particles are distributed uniformly throughout the entire simulation domain, and their velocities are drawn from the same Maxwell–Boltzmann velocity distribution as the main particles so as to ensure that the transport coefficients of both are the same [43]. However, bath particles do not participate in reactions and their sole purpose is to generate a constant desired diffusion of the main particles throughout the entire simulation domain.

N

vm

(2)

m=1

where Nξ is the number of particles in the cell and vm denotes the velocity of the mth particle. In order to calculate the post-collision velocity of the ith particle in subvolume ξ the following equation is used

vi (t + ) = V + ^ (vi (t )

(3)

V)

τ is the length of the time interval. In the all-species collision rule, the species of the particles are not distinguished, and the collision cell ξ is not the same subvolume for every time step as a result of the grid shifting. For a system that contains different chemical species that have different diffusion coefficients, the collision step usually involves both an all-species rule as well as a single-species rule. In this case we use a subscript for each of the different species α = 1,2,..,Nsp, where Nsp is the number of different species in the system, and the net effect of both the all-species and single-species rules can be written as [45]:

vi (t + ) = V + ^ (V

V ) + ^ ^ (v i (t )

V )

(4)

Vξ denotes the pre-collision average velocity of all particles in the subvolume ξ, while V is the average velocity of all particles of species α in the subvolume. The random rotation operators ^ and ^ act on all particles or all particles of species α in the subvolume, respectively. In the RMPC, a reaction Rj can be represented in the following general form: Nsp

j

kj

X

Nsp

=1

¯ jX

3. Inhomogeneous Stochastic Simulation Algorithm The ISSA is the augmentation of the Gillespie algorithm to generate exact realizations of the RDME [3,4,47]. The idea of the RDME is to divide a system, having Nsp species and M reactions, into NV cubic uniform subvolumes. Each of these subvolumes is assumed to be a wellstirred system of edge length h. Each of the M reactions Rj (j = 1, 2,..,M) may occur in every subvolume. As a result, the total number of reactions in the system amounts to M × NV in a way that Rjk represents reaction Rj in subvolume Ωk. ajk(x*k) denotes the propensity of Rjk where xik is the number of molecules of species Si in subvolume Ωk in the current state of the system. When reaction Rjk fires, the state of the system changes from x to x+vjk .ek. vjk is an Nsp dimensional column state change vector due to reaction Rjk and ek is an NV dimensional row vector whose entries are all zero except for the kth one which is unity. According to the ISSA, diffusion is assumed to be a unimolecular reaction, Rdikl, in which one molecule of type Si transfers between two d (xi *) represents the propensity neighboring subvolumes k and l. aikl function of this diffusive event. This propensity is equal to Di.xik /h2 with Di being the diffusion coefficient for species Si. When a unimd d olecular reaction Rikl fires, the system state x is updated to x + Ei ikl . Ei is defined as an Nsp dimensional column vector with its ith element d being 1 and all others are zero. ikl as an NV dimensional row vector d denotes the state change vector due to reaction Rikl . At each step of the ISSA, two random numbers r1 and r2 are computed from the uniform distribution in (0, 1). If a0 denotes the sum of the propensities for all the reactions and diffusive events in the system, then the time to the next event, δt, may be computed as

(5)

=1

where and ¯ are the stoichiometric coefficients for reactants and products for species α respectively. The kj is the rate at which reaction Rj takes place. The simulation algorithm treats the reactions as BirthDeath stochastic events [7] that take place at discrete time steps, just like collisions. The probability that reaction Rj takes place in a given subvolume ξ during a time interval of length τ is given by j

P j (N ) =

j

aj a

(1

e

)

a

(6)

The elements of the vector N = (N , N , …) consist of the number of particles in subvolume ξ for the species type as indicated by the superscript, and 1

a j = kj (V )

2

N j! (N j

j

)!

(7)

where the reaction rates kj(Vξ) are scaled by the volume Vξ of a subvolume, and a 0 = j a j [46]. Free-streaming updates the position of all particles at time t + τ using the post-collision and post-reaction velocities of a given particle according to: (8)

ri (t + ) = ri (t ) + vi (t + )

The transport properties of a system can be calculated using the RMPC dynamics. Green–Kubo relations in a 3 dimensional simulation [21] lead to a diffusion coefficient in terms of RMPC parameters as follows:

K T D= B m 2 (

3 1)(1

cos )

1

t=

1 1 ln a 0 (x) r1

(10)

with a0 being M

(9)

j=1 k=1

with ρ representing the average number of particles in a subvolume or collision cell, and ψ being the collision angle for the all-species collision rule. Eq. (9) however is applicable only if the following assumptions apply: first, variations in number density in the subvolume are negligible which is known to be valid when ρ ≥ 5, and second, the dimensionless mean free path = (KB T / m) / h must be greater than 0.6. Here, h is the length of an assumed cubic subvolume. There are two issues regarding the use of RMPC to simulate reaction

Nsp NV NV

NV

a0 (x) =

ajk (x*k) +

i =1 l=1 k=1

d aikl (xi *)

(11)

The next event is predicted as follows: if the index (j, k) is the smallest such that j

k

aj k j =1 k =1

(x*k ) > r2 a0 (x)

Then the reaction Rjk occurs first. If, however, the index obeys 25

(12)

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

(Intel 3.60 GHz quad-core CPU with 10GB RAM). 5.1. Diffusion The accuracy of the hybrid approach is first tested on the diffusion of particles in a computational domain [0, 1]. The constant diffusion coefficient is taken to be D = 1 mm2/s, with 2000 particles placed uniformly in the left half of the domain. Fully reflective boundary conditions are applied at x = 0 and x = 1. The first half of the domain is allocated to the ISSA and the second half is simulated by the RMPC; i.e. the interface is located at x = 0.5. In order to exploit the efficiency of the ISSA, the first half of the domain is divided into NISSA = 10 cubic cells of size hISSA = 0.05 mm. The RMPC domain is divided into Nx = 50 cubic subcells of size h = 0.01 mm. The number of subcells in the y and z directions are Ny = 5 and Nz = 5, respectively. The discrete time step τ for the RMPC simulation is set at 10−4 s from which the dimensionless mean free path is calculated to be λ = 1.414. In order for Eq. (9) to be valid, bath particles are placed uniformly throughout the RMPC domain so as to maintain a minimum bath density (ρbath = 5). The corresponding partial differential equation for the normalized density Z of molecules is given by

Fig. 1. Schematic illustration of the Ghost Cell Method. M

NV

aj j =1 k =1

k

(x*k ) +

i

l

k

i =1 l =1 k =1

The diffusion event

Rdikl

( )

aidk l x i * > r2 a 0 (x)

(13)

fires first.

4. The Ghost Cell Method for the RMPC–ISSA hybrid method In the Ghost Cell Method, the particles are allowed to cross over the interface in accordance with the compartment-based approach i.e. the ISSA. Particles are placed throughout the entire computational domain as prescribed by a chosen initial condition for the system. The time until the next (mesoscopic) ISSA event is compared to the fixed time step size of the (microscopic) RMPC algorithm, the earlier of which will take place. The particles in the ghost cells, which define a layer of cells in the microscopic subdomain adjacent to the interface, are prescribed by the microscopic (RMPC) dynamics. However, during the mescoscopic event, particles are allowed to jump over the interface from the mesoscopic (ISSA) subdomain into the ghost cells and vice versa in accordance to the ISSA. The details of the GCM is schematically illustrated in Fig. 1. The GCM has been implemented for the RMPC–ISSA hybrid method using the algorithm in Table 1. Governing part of the microscopic domain by the mesoscopic algorithm is justified by the fact that both algorithms are modeling the same diffusion phenomenon [28].

2Z Z =D 2 t x

0
(14)

which is solved subject to the initial condition Z(x,0) = 2H(0.5 − x) and boundary conditions Zx(0,t) = Zx(0,1) = 0. Here, H represents the standard Heaviside function. Applying separation of variables with D = 1 mm2/s, leads to the series solution

Z(x , t ) = 1 + n=1

4 n sin n 2

exp( n2 2t )cos(n x )

(15)

The analytical solution is compared to the hybrid model in Fig. 2, where very good agreement is seen. Note that the minimum number of particles in the RMPC cells was used. In order to suppress the stochastic noise of the simulation, the average of 100 realizations of the hybrid model is compared with the full RMPC and the analytical solution in Fig. 3. It should be noted that when modeling the whole domain using the RMPC, the number of bath particles initially placed in each cell needs to be increased. The result shown in Fig. 3 has been obtained with ρbath = 10. It can be seen in Table 2 that for the given discretization, the hybrid algorithm is more than 3 times faster than the full RMPC. As mentioned previously, mass, momentum and energy are

5. Numerical results In this section, the RMC-ISSA Hybrid algorithm is applied to three reaction-diffusion systems. In all cases, hybrid simulation results are compared to the solution of the corresponding partial-differential equations. Hybrid simulations were carried out on a standard desktop Table 1 Algorithm for RMPC–ISSA hybrid model. 1- Discretize the computational domain into subvolumes. 2- Specify the meso and micro subdomains and the interface separating them 3- Initialize the ISSA by placing particles in the meso subdomain in accordance with the given initial condition 4- Initialize the RMPC dynamics by placing particles in the micro subdomain in accordance with the given initial condition and choosing the discrete time step τ 5- Uniformly distribute the required number of bath particles throughout the micro subdomain 6- Initialize the time t:=0, tMPC=τ 7- Calculate the time for the next event in the meso subdomain and the ghost cells; δt using Eq. (10) 8- If t+δt < tMPC, change the state of the system in the meso subdomain or the ghost cells using the ISSA. t:=t+ δt Else Initialize the particles that have crossed over the interface into the ghost cells in preceding ISSA events by assigning them random position and velocity as per RMPC dynamics. Change the state of the system in the micro subdomain using the RMPC considering the fact that particles reflect upon hitting the interface. Count and store the number of the particles in the ghost cells for use in the next ISSA event. t=tMPC, tMPC=tMPC+ τ 9- Repeat steps 7 and 8 until the desired end of the simulation.

Fig. 2. The hybrid model simulation of a system where 2000 particles were initially distributed uniformly in the first half of the domain (t = 0.1 s). 26

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

Fig. 3. The average of 100 realizations of hybrid model simulations of a system where 2000 particles were initially distributed uniformly in the first half of the domain.

5.2. Morphogen gradient

Table 2 Run time speed-ups of the reaction diffusion systems. Reaction-diffusion system

Speed-up over the full RMPC algorithm

Simple diffusion Morphogen Gradient Schnakenberg system

3.25 6.15 35.4

The second example is the simulation of a diffusing Morphogen [51]. The computational domain is assumed to be semi-infinite [0, ∞). There is a source at x = 0 which produces molecules at a constant rate of kp = 6.25 s−1. The molecules are allowed to diffuse with D = 0.5 mm2/s. There is also a unimolecular depletion reaction at the rate of k = 1 s−1 that occurs through the whole computational domain. The governing differential equation for such a system can be expressed as follows: n (x , t ) t

D

=D

n (x , t ) x

=

2n (x , t )

x2

kp

ka (x , t ) x=0

x>0 (16)

with n(x,t) denoting the concentration of the particles. Due to the existence of a production source at x = 0, the number of molecules tends to decrease as they diffuse away from the source. For this reason, to utilize the efficiency of the compartment based modeling, a portion of the domain close to the source with a high concentration of molecules is simulated using the ISSA, while the rest of the domain having a low concentration of particles, is modeled using the RMPC where each individual particle can be tracked. The interval [0, 3] is considered to be the simulation domain. The interface separating the two regimes is placed at x = 2.5. The ISSA subdomain is discretized into NISSA = 250 cubic cells of size hISSA = 0.01 mm. The RMPC subdomain is discretized into Nx = 50, Ny = 1 and Nz = 1 cubic subcells of size h = 0.01 mm in the x, y and z directions, respectively. The time step for the RMPC is chosen to be 10−4 s resulting in a dimensionless mean free path of λ = 1. Fig. 5 shows the density of the particles at t = 3 s compared to the PDE solution, full RMPC and ISSA algorithms. As can be seen, all different approaches are in good agreement with one another. Given that the ISSA is the fastest of all the three numerical methods, this result may question the usefulness of the Hybrid method over the ISSA. In order to address this issue, the Morphogen gradient problem is solved once more with kp = 1.875 s−1 i.e. with a lower molecule production resulting in a system with a lower number of particles. Fig. 6a draws a comparison between the three methods in which the ISSA is seen to have failed. An enlarged view of the results close to the production source (x = 0) has been shown in Fig. 6b. Considering the fact that in

Fig. 4. The average of 100 realizations of hybrid model simulations of the system comparing slip versus no-slip BCs.

conserved in the collision step of the RMPC algorithm, however not all reactions are conservative. As such, MPC (RMPC without reactions) can lead to proper hydrodynamics unlike other particle-based algorithms. Additionally, RMPC can subsequently incorporate both slip and no-slip boundary conditions which cannot be implemented with ISSA. The need for slip versus no-slip for diffusive systems is discussed in [48–50]. Fig. 4 compares slip versus no-slip BCs in the RMPC subdomain for in a diffusion system. 27

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

reality the problem domain is semi-infinite, this failure is due to truncating the domain that for a system of low concentration introduces errors associated with small concentrations at large values of x. With the current discretization, as shown in Table 2, the Hybrid algorithm is seen to have sped up the simulation more than six fold. It needs to be noted that the purpose of modeling the above-given systems i.e. diffusion and Morphogen Gradient is just to test the validity as well as advantage of the proposed hybrid algorithm over the ISSA. The speed-up of the algorithm compared to the full RMPC, as will be seen in the next example, can be significantly increased through proper discretization of the domain. 5.3. Schnakenberg system Our third and last example simulated the diffusion of two reacting chemical species u and v with diffusion constants Du = 1 mm2/s and Dv = 0.1 mm2/s, respectively, in a computational domain of [0,1]. The left boundary (x = 0) is assumed to be partially adsorbing, while x = 1 is treated as fully reflective. The reactive mechanism follows the Schnakenberg kinetics [52,53]:

Fig. 5. Comparison between ISSA, Hybrid and full RMPC algorithms for a high molecule production rate for an average of 1000 realizations (t = 3 s).

Fig. 6. (a) Comparison between ISSA, Hybrid and full RMPC algorithms for a low molecule production rate for an average of 1000 realizations (t = 3 s); (b) An enlarged view of the comparison between ISSA, Hybrid and full RMPC algorithms for a low molecule production rate.

28

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

Fig. 7. The hybrid simulation of a reaction-diffusion Schnakenberg system for species u (t = 0.1 s); (b) The hybrid simulation of a reaction-diffusion Schnakenberg system for species v (t = 0.1 s). c1

0 U 0

c2 c3

0, c2 = 6 s

(18)

1

2U + V

3U , c4 = 3 mm2/(mol2.s)

(20)

The Schnakenberg model is capable of pattern formation, and U is often referred to as the activator, while V is the inhibitor. The corresponding reaction-diffusion system is given by 2u u = Du 2 + c1 t x

c2 u + c4

2v v = Dv 2 + c3 t x

c4 u2v

u2v

Pa =

v (x , 0) = 1

2

3.3321 K

(24)

t 2D

2

+ 3.35669

3

1.52092

4

(25)

with K describing the reactivity of the boundary. The where = same probability Pa can be used with the RMPC provided that the dimensionless mean free path of the particle in the RMPC, ( (KB T /m ) ), is set equal to the root mean square displacement of the particle in Brownian Dynamics ( 2D ). The details of modeling reactive boundary conditions using the RMPC can be found in reference [43]. It should also be noted that for such a system with more than one species, the RMPC involves only a single species collision for each species with a different value of KBT/m . The interface is placed at x = 0.2. This is

(21) (22)

Subject to the following initial and boundary conditions:

u (x , 0) = 1,

v = 10v (0, t ) x

In a particle-based method, one can simulate a partially adsorbing boundary condition at x = 0 as follows: If the position of the particle becomes negative, then the particle is adsorbed (or rather removed from the system) with probability Pa, and reflected otherwise. Using Brownian dynamics, Andrews [54] provides an estimate for Pa as follows:

(19)

V , c3 = 8 mol/(mm.s) c4

u = 10u (0, t ), x

(17)

U , c1 = 2 mol/mm.s

(23) 29

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

Fig. 8. The average of 100 realizations of a reaction-diffusion Schnakenberg system for species u; The average of 100 realizations of a reaction-diffusion Schnakenberg system for species v.

because only a small part of the domain close to the partially reactive boundary i.e. x = 0, where individual particles need to be tracked, is modeled using the RMPC, whereas the rest of the domain is governed by the ISSA. The RMPC subdomain is divided into Nx = 20 cubic subcells of size h = 0.005 mm. The number of cells in the y and z directions are Ny = 4 and Nz = 4, respectively. The ISSA subdomain is discretized into 45 cubic cells of size hISSA = 0.02 mm. Initially, 12,800 particles are uniformly distributed throughout the entire domain. This uniform distribution corresponds to 256 particles in each ISSA cell and 4 particles in each RMPC subcell. The RMPC time step is set to be 10−4 s leading to dimensionless mean free paths of λu = 2.828 and λv = 0.894 for species u and v, respectively. Figs. 7a and b show results compared to the PDE solution for species u and v, respectively at t = 0.1 s. The average of 100 realizations of the hybrid algorithm has been compared with the full RMPC algorithm in Figs. 8a and b for the two different species. The use of the hybrid algorithm as reported in Table 2, has resulted in a speed-up of more than 35 times over the full RMPC

algorithm. The speed-ups listed in Table 2 strongly depend on the very system being simulated, the way the domain is split between the RMPC algorithm and the ISSA and, more importantly, on the way each subdomain is discretized. One can always increase the speed-up by either increasing the size of the ISSA subdomain or coarsening the ISSA mesh. 6. Discussions The development of hybrid algorithms stems from the need for saving both computational time and accuracy when fine and microscopic details are required only in part of the problem domain. One class of such hybrid methods apply to micro–meso models. Unlike the previous micro–meso hybrid methods where the microscopic subdomain is simulated using either Brownian Dynamics or Green's Function Reaction Dynamics, the present work utilizes the less detailed and more efficient RMPC dynamics for modeling the trajectories of 30

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al.

individual particles and reactions. The mesoscopic subdomain is simulated using the ISSA. The interface between the two domains is modeled using the Ghost Cell Method. In order to maintain a fixed diffusion coefficient in the RMPC subdomain in the course of the simulation, bath particles are initially uniformly distributed throughout the whole simulation domain. These bath particles do not participate in reactions and their sole purpose is to generate a constant desired diffusion of the main particles throughout the entire simulation domain. The main advantage of RMPC in comparison with other methods such as Brownian Dynamics is that it can model hydrodynamic interactions due to conservation of mass, momentum and energy of the system in its free stream and collision steps and also the generation of fewer random numbers. The proposed hybrid method is tested on three benchmark problems. The first system is a simple diffusion with the interface placed at the middle of the domain. The second model is that of a Morphogen in which the computational domain away from the molecule production source, where there is a low copy number, is modeled with the microscopic algorithm. The third problem is the simulation of the Schnakenberg system in the presence of a reactive boundary, in which a small portion of the domain close to the reactive boundary is modeled with microscopic details. Unlike the RMPC subdomains that are smaller than the ISSA subcells, the ISSA subdomain has 1 cell in each of the y and z direction. This ensures greater accuracy in the RMPC subdomain, in addition to increasing the speed of the ISSA. Additionally, the boundaries in the y and z directions are treated as fully reflective in the RMPC subdomain. In all the problems studied, reasonable agreement with the deterministic solutions of the corresponding differential equation is achieved at reduced computational cost compared to the full RMPC algorithm. The discrepancies in results stem mainly from RMPC for which the simulated diffusion may not be the same as predicted by Eq. (9) as a result for the mean free path not being large enough. However, we showed that adding enough bath particles can provide much better agreement over a range of mean free paths. Our results can be applied so long as the particles in the system are small enough so that volume exclusion effects are negligible and macromolecular crowding can also be neglected. Both the ISSA and the RMPC would need modifications if they are to incorporate these effects [55–57]. It is worth noting that generally the stochastic mean behavior of reactiondiffusion systems is comparable to the PDE solution for reactive systems governed by zeroth and first order reactions for which detailed balance holds [58]. In the Schnakenberg model used in this paper, the number of molecules present in the system has thus been chosen to be large enough so that a comparison to the PDE solution was possible.

[3] A.B. Sturzia, C.J. Lumsden, Stochastic simulation of coupled reaction–diffusion processes, J. Comput. Phys. 127 (1996) 196–207. [4] S.A. Isaacson, C.S. Peskin, Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations, SIAM J. Sci. Comput. 28 (2006) 47–74. [5] D. Fange, J. Elf, Noise-induced min phenotypes in E. coli, PLOS Comput. Biol. 2 (2006) 637–648. [6] S. Smith, R. Grima, Spatial stochastic intracellular kinetics: a review of modelling approaches, Bull. Math. Biol. (2018) 1–50. [7] K. Rohlf, S. Fraser, R. Kapral, Reactive multiparticle collision dynamics, Comput. Phys. Commun. 179 (2008) 132–139. [8] J. Elf, M. Ehrenbeg, Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases, Syst. Biol. 1 (2004) 230–236. [9] L. Ferm, A. Hellander, P. Lotstedt, An adaptive algorithm for simulation of stochastic reaction–diffusion processes, J. Comput. Phys. 229 (2010) 343–360. [10] J.M.A. Padgett, S. Ilie, An adaptive tau-leaping method for stochastic simulations of reaction-diffusion systems, AIP Adv. 6 (2016) 035217. [11] R. Strehl, S. Ilie, Hybrid stochastic simulation of reaction-diffusion systems with slow and fast dynamics, J. Chem. Phys. 143 (2015) 234108. [12] M. Barrio, K. Burrage, A. Leier, T. Tian, Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation, PLoS Comput. Biol. 2 (2006) 1017–1030. [13] A. Sayyidmousavi, S. Ilie, An efficient hybrid method for stochastic reaction-diffusion biochemical systems with delay, AIP Adv. 7 (2017) 125305. [14] R. Holley, The motion of a heavy particle in an infinite one dimensional gas of hard spheres, Probab. Theory Relat. Fields 17 (1971) 181–219. [15] R. Strehl, K. Rohlf, Multiparticle collision dynamics for diffusion-influenced signaling pathways, Phys. Biol. 13 (2016) 046004. [16] S.S. Andrews, D. Bray, Stochastic simulation of chemical reactions with spatial resolution and single molecule detail, Phys. Biol. 1 (2004) 137–151. [17] J.S. van Zon, P.R. ten Wolde, Greens-function reaction dynamics: a particle-based approach for simulating biochemical networks in time and space, J. Chem. Phys. 123 (2005) 234910. [18] J.S. van Zon, P.R. ten Wolde, Simulating biochemical networks at the particle level and in time and space: green's function reaction dynamics, Phys. Rev. Lett. 94 (2005) 128103. [19] A. Malevanets, R. Kapral, Mesoscopic model for solvent dynamics, J Chem Phys 110 (1999) 8605–8613. [20] A. Malevanets, R. Kapral, Solute molecular dynamics in a mesoscale solvent, J. Chem. Phys. 112 (2000) 7260–7269. [21] K. Tucci, R. Kapral, Mesoscopic model for diffusion influenced reaction dynamics, J. Chem. Phys. 120 (2004) 8262–8270. [22] K. Tucci, R. Kapral, Mesoscopic multiparticle collision dynamics of reaction-diffusion fronts, J. Phys. Chem. B 109 (2005) 21300–21304. [23] C.A. Smith, C.A. Yates, Spatially extended hybrid methods: a review, J. R. Soc. Interface 15 (2018) 20170931. [24] M.B. Flegg, S.J. Chapman, R. Erban, The two regime method for optimizing stochastic reaction–diffusion simulations, J. R. Soc. Interface 9 (2012) 859–868. [25] M. Robinson, M.B. Flegg, R. Erban, Adaptive two-regime method: application to front propagation, J. Chem. Phys. 140 (2014) 124109. [26] M.B. Flegg, S.J. Chapman, L. Zheng, R. Erban, Analysis of the two-regime method on square meshes, SIAM J. Sci. Comput. 36 (2014) 561–588. [27] A. Hellander, S. Hellander, P. Lotstedt, Coupled mesoscopic and microscopic simulation of stochastic reaction–diffusion processes in mixed dimensions, Multiscale Model. Simul. 10 (2012) 585–611. [28] M.B. Flegg, S. Hellander, R. Erban, Convergence of methods for coupling of microscopic and mesoscopic reaction–diffusion simulations, J. Chem. Phys. 289 (2015) 1–17. [29] J.M. Yeomans, Mesoscale simulations: lattice Boltzmann and particle algorithms, Physica A 369 (2006) 159. [30] P.J. Hoogerbrugge, J.M.V.A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, EPL 19 (1992) 155. [31] T. Ihle, D.M. Kroll, Stochastic rotation dynamics: a Galilean-invariant mesoscopic model for fluid flow, Phys. Rev. E. 63 (2001) 020201. [32] T. Ihle, D.M. Kroll, Stochastic rotation dynamics: I. Formalism, Galilean invariance, and Green–Kubo relations, Phys. Rev. E 67 (2003) 066705. [33] T. Ihle, D.M. Kroll, Stochastic rotation dynamics: II. Transport coefficients, numerics, and long-time tails, Phys. Rev. E 67 (2003) 066706. [34] H. Noguchi, G. Gompper, Transport coefficients of off-lattice mesoscale-hydrodynamics simulation techniques, Phys. Rev. E 78 (2008) 016706. [35] E. Allahyarov, G. Gompper, Mesoscopic solvent simulations: multiparticle-collision dynamics of three-dimensional flows, Phys. Rev. E 66 (2002) 036702(9). [36] K. Mussawisade, M. Ripoll, R.G. Winkler, G. Gompper, Dynamics of polymers in a particle-based mesoscopic solvent, J. Chem. Phys. 123 (2005) 144905(11). [37] N. Kikuchi, J.F. Ryder, C.M. Pooley, J.M. Yeomans, Kinetics of the polymer collapse transition: the role of hydrodynamics, Phys. Rev. E 71 (2005) 061804(8). [38] J.T. Padding, A.A. Louis, Hydrodynamic interactions and Brownian forces in colloidal suspensions: coarse-graining over time and length scales, Phys. Rev. E 74 (2006) 031402(29). [39] D.J. Earl, C.M. Pooley, J.F. Ryder, I. Bredberg, J.M. Yeomans, Modeling microscopic swimmers at low Reynolds number, J. Chem. Phys. 126 (2007) 064703(10). [40] G. Rückner, R. Kapral, Chemically powered nanodimers, Phys. Rev. Lett. 98 (2007) 150603(4). [41] K. Rohlf, Stochastic phase-space description for reactions that change particle numbers, J. Math. Chem. 45 (2009) 141–160. [42] J. Lipková, K.C. Zygalakis, S.J. Chapman, R. Erban, Analysis of Brownian dynamics simulations of reversible bimolecular reactions, SIAM J. Appl. Math. 71 (2011)

Acknowledgments The authors wish to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for the support of this research. Special thanks also go to the reviewers for their comments which resulted in substantial improvement of the manuscript. Funding This work was supported by funds from the Ryerson University Mathematics Department and the Faculty of Science, as well as a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery grant (Grant number: RGPIN-2017-04672, RGPIN-201505723). References [1] D.T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem. 81 (1977) 2340–2361. [2] D.T. Gillespie, Approximate accelerated stochastic simulation of chemically reacting systems, J. Chem. Phys. Chem. 155 (2001) 1716–1733.

31

Mathematical Biosciences 312 (2019) 23–32

A. Sayyidmousavi, et al. 714–730. [43] A. Sayyidmousavi, K. Rohlf, Reactive multi-particle collision dynamics with reactive boundary conditions, Phys. Biol. 15 (2018) 046007. [44] H. Hijar, G. Sutmann, Hydrodynamic fluctuations in thermostatted multiparticle collision dynamics, Phys. Rev. E 83 (2011) 046708. [45] R. Kapral, Multiparticle Collision Dynamics: simulation of complex systems on mesoscales, Adv. Chem. Phys. 140 (2008) 89–146. [46] J.X. Chen, J.X. Zhu, Y.H. Zhau, W.G. Sun, J.R. Xu, H.P. Ying, Simulating bistable biochemical systems by means of reactive multiparticle collision dynamics, Commun. Nonlinear Sci. Numer. Simul. 19 (7) (2014) 2505–2512. [47] I. Hepburn, W. Chen, S. Wils, E. De Schutter, Steps: efficient simulation of stochastic reaction–diffusion models in realistic morphologies, BMC Syst. Biol. 6 (2012) 36. [48] V. Mortazavi, M. Nosonovsky, Friction-induced pattern formation and Turing systems, Langmuir 27 (2011) 4772–4779. [49] D.A. Garzon-Alvarado, C.H. Galeano, J.M. Mantilla, Computational examples of reaction-convection-diffusion equations solution under the influence of fluid flow: a first example, Appl. Math. Model. 36 (2012) 5029–5045. [50] J. Wei, M. Winter, Flow-distributed spikes for Schnakenberg kinetics, J. Math. Biol.

64 (2012) 211–254. [51] F. Tostevin, P.R. ten Wolde, M. Howard, Fundamental limits to position determination by concentration gradients, PLoS Comput. Biol. 3 (2007) 0763–0771. [52] J. Schnakenberg, Simple chemical reaction systems with limit cycle behavior, J. Theor. Biol. 81 (1979) 389–400. [53] R. Erban, S.J. Chapman, Reactive boundary conditions for stochastic simulations of reaction-diffusion processes, Phys. Biol. 4 (2007) 16–28. [54] S.S. Andrews, Accurate particle-based simulation of adsorption, desorption, and partial transmission, Phys. Biol. 6 (2009) 046015. [55] C. Cianci, S. Smith, R. Grima, Capturing Brownian dynamics with an on-lattice model of hard-sphere diffusion, Phys. Rev. E 95 (2017) 052118. [56] S. Smith, C. Cianci, R. Grima, Macromolecular crowding directs the motion of small molecules inside cells, J. R. Soc. Interface 14 (131) (2017) 20170047. [57] C. Echeveria, K. Tucci, R. Kapral, Diffusion and reaction in crowded environments, J. Phys. 19 (2007) 065146–065158. [58] S. Smith, C. Cianci, R. Grima, Analytical approximations for spatial stochastic gene expression in single cells and tissues, J. R. Soc. Interface 13 (2016) 20151051.

32