Stochastic simulations of the Schnakenberg model with spatial inhomogeneities using reactive multiparticle collision dynamics

1.   Introduction

Systems in which reactions and diffusion occur together have been of interest for many decades. It has long been known that diffusion can alter the reactive dynamics in the absence of diffusion, leading to diffusive instabilities such as Turing patterns. Studying the effects that noise may have on the system dynamics has more recently led to a wealth of stochastic simulation studies. Stochastic simulation methods are also necessary whenever the ODE or PDE approximation ceases to hold, which can happen whenever there are low concentrations of a particular chemical species limiting when reactions involving that chemical species can actually take place locally.

Several stochastic simulation algorithms exist, the most popular of which are based on the stochastic simulation algorithm (SSA) initially introduced by Gillespie ^[1,2]. SSA has more recently been extended to allow for incorporation of diffusion of the chemical species, as well as spatial inhomogeneities in what is referred to as the inhomogeneous stochastic simulation algorithm (ISSA) ^[3,4,5]. There are many benefits of the ISSA method, in particular the computational efficiency when the chemical species have large concentrations. Shortcomings for regions where there may be low concentrations of a particular chemical species are addressed through various extensions of the ISSA using different approaches. These include adaptive time stepping methods ^[6,7,8], or hybrid approaches ^[9,10,11,12]. Several key drawbacks of these models are the inability to consider velocity-dependent reaction rates which may be more realistic when discussing reactive dynamics in molecular self-assembly for example, or for incorporating slip versus no-slip boundary conditions when simulating realistic biochemical networks for which this may be important. Additionally, individual particle-tracking is also not possible with traditional SSA or ISSA models, nor their direct extensions that include the next subvolume method ^[13].

Microscopic models that incorporate stochastic dynamics comparable to ISSA with more details at the particle-level include Brownian dynamics (BD) simulations. Brownian dynamics (BD) models provide more efficient dynamics at a more detailed level compared to ISSA. In BD, particle trajectories can be traced, however velocities of individual particles are not modeled explicitly. Recent popular simulations tools for BD methods include Smoldyn ^[14] and Green's function reaction dynamics ^[15,16]. Extensions to incorporate particle velocities into the dynamics has recently been used in the context of reactive boundary conditions ^[17], giving rise to so-called Langevin dynamics for reaction-diffusion processes. A nice overview of velocity jump approaches applied to reaction-diffusion applications in biological systems can be found in ^[18].

Turing patterns for the Schnakenberg model have been produced in stochastic models using the traditional RDME, as well as a modified RDME approach ^[19]. An example of a Turing pattern was provided that formed using a modified RDME approach, where the two chemical species evolve on different mesh sizes. The traditional RDME simulation did not lead to a Turing pattern nor the comparable PDE solution, while the modified RDME approach did. In ^[20], Turing patterns were found using a generalized Smoluchowski framework, and compared to a PDE solution using a perturbed heterogeneous form for one of the reaction rates so as to model more accurately the boundary effect of the diffusion-controlled reactions of the Smoluchowski model.

More recent particle-based methods that are capable of incorporating stochastic effects for reaction-diffusion systems include the extension of the multi-particle collision (MPC) dynamics for particle-based fluid flow or diffusive motion, to include reactions, giving rise to reactive multi-particle collision (RMPC) dynamics. The MPC dynamics was first introduced in ^[21,22], and extended to include reaction mechanisms thereafter ^{[23,24,25,26]}. MPC/RMPC is a more efficient simulation strategy than comparable Lattice-Boltzmann ^[27] or dissipative particle dynamics ^[28], in that particle positions and velocities are not constrained to a lattice, but rather are three-dimensional vectors with components in $\Re$. Incorporation of velocities in addition to positions allows for the possibility to model diffusive as well as fluid behaviour. Most MPC publications exploit the ability of the stochastic particle-based method to simulate fluid flow ^[29,30,31]. Much work has been done in the development of connecting transport coefficients to averages of stochastic simulation parameters, including theoretical expressions for the diffusion coefficient ^{[21,22,32,33,34,35,36,37]}. However, more recent work focuses on developing RMPC for stochastic simulations of biochemical networks that are more traditionally modeled using ISSA or BD methods ^{[25,26,37,38]}. Limitations and capabilities of RMPC dynamics applied to exploit the diffusive nature in the reactive system is less understood and less developed, and this paper focuses on exploring and highlighting the diffusive aspect of RMPC in a reaction-diffusion system, thus setting the stage for future work that can incorporate slip versus no-slip boundary conditions in finite domains, or velocity-dependent reaction rates that might be important for some applications ^[39,40,41]. To the best of our knowledge, RMPC simulations have not yet been reported in the literature to have generated Turing patterns.

The paper is organized as follows. In Section 2 the RMPC dynamics is presented, as well as how the parameters of the simulation method can be fine-tuned to specify a constant diffusion coefficient on the particle system that is comparable to Brownian dynamics. The concept of how bath particles can be effectively introduced into the system dynamics to fix the diffusion coefficient value when there may be large spatial inhomogeneities in the system is also introduced in this section. Two different ways to vary the diffusion coefficient are used (ⅰ) varying $k_B T/m$, and (ⅱ) free-streaming with a certain probability and varying the probability to free-stream. Numerical results for both the diffusion-only case (eg. no reactions), and for the full reaction-diffusion system are presented in Section 3. Solutions are compared to the corresponding PDE solution for a spatially inhomogeneous initial condition with and without bath, where the stochastic dynamics is significantly different without bath particles. Next, the concept of free-streaming with a certain probability is used in long-time simulations, and shown to produce Turing patterns. A summary of the results as well as important conclusions and future work are given in the final Section 4.

2.   Reactive multi-particle collision dynamics

The system contains $N$ particles in the simulation volume $V$ which is divided into $N_c$ cells each with length $dx$, $dy$ and $dz$ in the $x$, $y$, and $z$-directions respectively. Taking $L_x$, $L_y$ and $L_z$ cells in the directions indicated by the respective subscripts, it follows that $V = (L_x dx)(L_y dy)(L_z dz)$. At discrete time steps $\tau$, the particles undergo collisions and reactions, followed by free-streaming. Reactions and collisions occur separately and independently in each cell, and particles move ballistically in the free-streaming step.

For the collisions, velocities of the $N$ particles in the system are updated based on the center of mass velocity in the collision cell in which they reside at that time. All particles are assumed to have the same mass $m = 1$. Suppose $\xi = 1, \ldots N_c$ is the cell number label, and $i = 1, \ldots N$ the particle number. Suppose particle $i$ has velocity ${\bf{v}}_i$, is located at position ${\bf{x}}_i$ and is of species $\alpha_i$ for a system containing $\alpha = 1, \ldots N_{sp}$ different species. Then the single-species rule is applied to the particles according to

where ${{\bf{V}}_\xi'}^{\alpha_i'}$ is the center of mass velocity of all particles with the same species, and in the same cell $\xi$, as particle $i$, ${\bf{v}}_i'$ is the pre-collision velocity of particle $i$, and ${\bf{v}}_i$ is the velocity of the particle after the collision rule has been applied. The rotation operator ${\hat{\omega}_\xi}^{\alpha}$ is a rotation matrix in the form

with $\psi = \psi_\alpha$, and $l_x = \sqrt{1-\theta^2} \cos{\phi}$, $l_y = \sqrt{1-\theta^2} \sin{\phi}$, $l_z = \theta$, $\phi$ and $\theta$ are random numbers on the intervals $[0, 2\pi]$ and $[-1, 1]$, respectively. Here $s_\psi = \sin{\psi}$ and $c_\psi = \cos{\psi}$, and at each time step and for each cell $\xi$, $\psi$ is chosen randomly from the set $\Omega$.

In order to maintain a constant diffusion coefficient throughout the system volume, the rotation angle $\psi_\alpha$ is varied at each time step for each species $\alpha$ in the system, and calculated from

where $\rho_\alpha$ is the average number of particles of species $\alpha$ in a cell based on a system average. For this choice in $\psi_\alpha$, the diffusion coefficient for that species becomes ^{[32,33,34,35,36,38]}

which is then analogous to Brownian dynamics ^[38]. Note that equation (2.3) leads to admissible values for $\psi_\alpha$ if and only if the right hand side is in the interval $[-1, 1]$, which will be the case if and only if $\rho_\alpha > 4$.

The numerically efficient choice of $\rho_{\alpha}$ comes from a global system average, that is from the total number of particles of species type $\alpha$ in the system divided by the total number of collision cells. This method was first introduced in ^[38]. For this choice of rotation angle, it can be seen from (2.4) that the value for $k_B T/m$ determines the value of the diffusion coefficient for that species. As such, different species in the system have different diffusion coefficients as specified by their $k_B T/m$ value.

In this work, only the single-species collision rule is used, and momentum exchange between different species occurs through reactions only. This is expected to be a good approximation as long as the number of reactive events that lead to momentum exchange between different species is sufficiently small.

In order to minimize correlations arising from particles residing in any given cell for several time steps as a result of small velocities of particles, standard grid-shifting methods are applied.

In the event that there are large spatial variations in particle density, it may be necessary to add bath particles so as to avoid large spatial variations in local diffusivities of the particles. For example, if the local number density $\rho^\xi_\alpha$ is substantially different from the global average $\rho_\alpha$, there may be large enough differences between the $D_\alpha$ value from (2.4) as compared to the local diffusivity

the latter of which would manifest itself as spatial variations of the diffusion coefficient. In that case, bath particles for the species $\alpha$ are added uniformly throughout the system volume and partake in the collision and free-streaming steps only. Bath particles do not take part in any reactions, and their sole purpose is to ensure that the local number density $\rho_\alpha^\xi$ of a given species is as close as possible to its global average. Figure 1 shows the effect that changes in number density has on the local diffusion coefficient for a given collision angle. Additionally, at small mean-free paths $\lambda$, where $\lambda = \sqrt{\frac{k_B T}{m}} \frac{\tau}{h}$, grid-shifting may not be effective so that the molecular chaos approximation may no longer hold even with grid-shifting, giving rise to deviations from the theoretical expression (2.4) as a result of correlations ^[36]. Theoretical expressions for the diffusion coefficient in those cases do not exist, but a trial-and-error approach can be used to find the best bath so as to obtain the desired value for the diffusion coefficient with collision angle chosen from the provided formula.

Figure 1.  Plot of local diffusion coefficient $D_\alpha^\xi$ as a function of $\rho_\alpha^\xi$ for different collision angles $\psi_\alpha$ as calculated from a global approach. The limiting case for $\psi_\alpha^{max}$ corresponds to the case when $\cos{\psi_\alpha} = \lim\limits_{\rho_\alpha \rightarrow \infty} \frac{\rho_\alpha+2}{2(1-\rho_\alpha)} = -\frac{1}{2}$. For this figure, $k_B T/m = 2000 \ mm^2/sec^2$.

DownLoad: Full-Size Img PowerPoint

Glancing at Eq. (2.5), changing the diffusion coefficient is most naturally achieved by changing the value for $k_B T/m$, especially when dealing with more than one chemical species in the system. Simulations using this approach are carried out in Section 3.1.1 for the diffusion of a highly concentrated peak in the center of the simulation domain, and for the reaction-diffusion system subject to similar initial conditions for both chemical species in the system (Section 3.2.1). Another way one can achieve different diffusivities, is to free-stream with a given probability. That is, if a random number drawn is less than a specified probability $p$, then the particle free-streams, and if not, then its position remains unchanged for that time step. In this way, one can simulate different diffusivities $pD^*$ for a given choice in $k_B T/m$, through variations in $p$. Simulations demonstrating this technique are provided in Section 3.1.2 for the diffusion of a highly concentrated initial condition, and used to generate Turing patterns in Section 3.2.2. The probability to free-stream approach was motivated by compartment based models where, in a given time step, either a reaction or a diffusive event happens, but never both in the same time step. In RMPC, reactions and free-streaming/diffusion takes place at every time step. Through the probability to free-stream approach, the same $k_B T/m$ value can be used for all chemical species in the system, and diffusion can be varied to a large degree, in a comparable way to compartment-based models. This is particularly advantageous for applications that require long simulation times, such as what is typically required for Turing patterns to form.

3.   Numerical results

Our methods are applied to simulate a reaction-diffusion mechanism where the reactions are governed by the Schnakenberg model. The reactive mechanism is given by

where $k_1^*$ to $k_4^*$ are the reaction rates for the reactions as indicated.

Stochastic simulations in a 3-dimensional domain are performed with diffusion measured in one direction, non-homogeneous initial conditions, and different diffusion coefficients for the two species.

The corresponding macroscopic governing partial differential equation for concentrations $u^*$ and $v^*$ (in number of particles per unit volume) for comparison is given by

where $D_u^*$ and $D_v^*$ are the diffusion constants of species $U$ and $V$ respectively. This system of PDEs is solved numerically subject to the boundary conditions

where $x^*$ is measured in $mm$ and $t^*$ in $sec$.

In terms of number of particles $U$, $V$ per (dimensionless) unit volume, where $u^* = U/ {h^*}^3$, $v^* = V/{h^*}^3$, the PDE can be written as

subject to the conditions

where $k_1^* = k_1 /{h^*}^3$, $k_2^* = k_2 /{h^*}^3$, $k_3^* = k_3$ $k_4^* = k_4 {h^*}^6$.

The deterministic solution for $U (x^*, t^*), V (x^*, t^*)$ can be compared to the stochastic local number densities $\rho_u^\xi, \rho_v^\xi$ from the RMPC simulations. For RMPC simulations, a globally averaged number density is used for calculation of the collision angle at each time step. More specifically, $\rho_u (t) \ = \ \sum_\xi \rho_u^\xi(t)/N_c$ where $N_c$ is the number of collision cells, is used in (2.4) and likewise for $\rho_v(t)$. This approach should be feasible as long as the globally computed number density $\rho_\alpha$ in (2.4) is comparable to the local number density $\rho^\xi_\alpha$ in (2.5) throughout the simulation volume. In the event that it is not, bath particles are added to correct for this.

The Schnakenberg model has often been used as an example to demonstrate Turing pattern formation. Turing patterns, in this case, would be a spatially inhomogeneous steady-state created as a result of diffusion-induced instability of the spatially homogeneous steady state. Necessary and sufficient conditions for Turing patterns to form for the Schnakenberg model in a one-dimensional domain are well-known, and generally require a large enough length, or a large enough ratio of the diffusivity of $V$ relative to that of $U$. To generate patterns, we used the same $k_B T/m$ value for the two chemical species in the system, and the probability of free-streaming concept to lower the probability of one chemical species relative to the other. This was done so that there would be no long-term effects of energy transfer from one species to the other through reactions (e. g., $2U+V \rightarrow 3U$) which could possibly change the diffusivity of $U$ in the long-term, as Turing patterns generally require very long simulation times. Furthermore, the long-term efficacy of bath particles is difficult to ascertain since a spatially homogeneous state is generally produced during relatively short simulation times compared to Turing pattern formation simulations. The bath particles serve to provide good agreement with the PDE solution in the initial transient dynamics reported here, but its long-term effectiveness for the length of simulation times required for Turing patterns to form has not been fully explored and tested here. Using the probability of free-streaming concept provides a simpler alternative for long-term simulations without requiring extensive computational studies to try to measure the long-term efficacy of bath particles.

3.1.   Diffusion only case

In all simulations for the case when the chemical species simply diffuses, the number of cells in the $x$, $y$ and $z$ directions are fixed at $N_x = 200$, $N_y = 1$ and $N_z = 1$ respectively. The dimensionful time step is $\tau^* = 10^{-4} \ sec$, and the grid size is $h^* = 5 \mu m$ with cubic cells. The initial condition described in Eq. (3.5) is used, and fully reflective boundary conditions are imposed at all boundaries of the simulation volume. Grid-shifting is also used.

3.1.1.   Varying $k_B T/m$

In order to test the use of a globally averaged number density in determining the collision angle to simulate an a priori diffusivity, stochastic MPC simulations (RMPC without reactions) are performed. There is only one species in the system, and the number of particles will satisfy the standard diffusion equation. The spatially inhomogeneous initial condition of (3.5) is used for the species. Different choices in $k_B T/m$ are considered, and the MPC results are compared to the numerical solution of the diffusion equation. Figure 2 shows the comparison of the MPC results to the PDE solution for various choices in $D$ when no bath is used. The MPC time step is $\tau^* = 10^{-4} \ sec$, and solutions are presented at $t = 0.01 \ sec$. For the larger diffusion coefficients, this smaller time is necessary to capture a transient state as the system reaches the spatially homogeneous steady state rather quickly. It can be seen that for small diffusion coefficients, the MPC solution differs significantly from the PDE solution, while for large values of $D$, there is good agreement. Recall that $k_B T/m$ is varied so as to achieve different diffusion coefficients, and in doing so, the mean free paths vary as well. Since for our simulations $D = \frac{k_B T}{m} \frac{\tau}{2}$, the mean-free path can be computed using $\lambda = \sqrt{\frac{k_B T}{m}} \frac{\tau}{h} = \sqrt{2D \tau}/{h}$.

Figure 2.  Comparison of the MPC solution (RMPC without reactions/diffusion-only case) to the PDE solution at $t^* = 0.01 \ sec$ for different mean-free paths without bath particles. $D = 0.1, 1, 2$ and $5 \ mm^2/sec$. Stochastic results are averages over 100 realizations.

DownLoad: Full-Size Img PowerPoint

Bath particles can be added to the system, such that the MPC results agree well with the PDE solution. Using a trial and error approach, Table 1 shows the bath density used to obtain good agreement with the PDE solution. The effectiveness of this approach can be seen in Figure 3 where the best bath density given in Table 1 is used in the MPC simulations. It is clear from the figure that very good agreement can be achieved for a large range of diffusion coefficients, and is in fact necessary to reach good agreement for the smaller diffusion coefficients considered.

Table 1.  Bath effectiveness as a function of mean free path $\lambda = \sqrt{2D \tau}/h$.

$\frac{k_B T}{m} = 2 \frac{D^*}{\tau^*}$	$\lambda$	$D^*$	Best bath, $\rho^{bath}$	No bath/Agreement
($mm^2/sec^2$)		($mm^2/sec$)
2000	0.894	0.1	220	very poor
20000	2.82	1	22	poor
40000	4	2	42	poor
100000	6.32	5	7	relatively good
2 x $10^5$	8.9443	10	7	relatively good

---

 | Show Table

DownLoad: CSV

Figure 3.  Comparison of the MPC solution (RMPC without reactions/diffusion-only case) to the PDE solution at $t^* = 0.01 \ sec$ for different mean-free paths with bath particles as per Table 2. $D = 0.1, 1, 2$ and $5 \ mm^2/sec$. Stochastic results are averages of 100 realizations.

DownLoad: Full-Size Img PowerPoint

3.1.2.   Free-streaming with probability $p$

Using a similar trial-and-error approach to find the best bath for comparison with the PDE solution, particles diffuse with probability $p$ for a fixed $k_B T/m$ value. The largest $k_B T/m$ value from Table 1 is chosen, and probabilities $p$ are assigned so that $p D^*$ provides the same range of diffusivities studied in the previous section. The simulation parameters used, and the required best bath are reported in Table 2. We found that lower bath particles were required for the simulations to give good agreement with the PDE solution, although no bath seemed able to lead to good agreement with the PDE solution in the smallest diffusivity ($pD^* = 0.1 mm^2/sec$). The number of particles as a function of $x$ for one realization for each of the Table 2 cases is plotted in Figures 4–6.

Table 2.  Varying Diffusion coefficient for a fixed $k_B T/m$ using probability $p$ to free-stream.

$\frac{k_B T}{m}$	$p$	$p D^*$	$\rho^{bath}$
($mm^2/sec^2$)		($mm^2/sec$)
2x$10^5$	1	10	7
2x$10^5$	0.5	5	14
2x$10^5$	0.2	2	16
2x$10^5$	0.1	1	10
2x$10^5$	0.01	0.1	5,200

---

 | Show Table

DownLoad: CSV

Figure 4.  Comparison of one realization of the MPC solution (RMPC without reactions/diffusion-only case) to the PDE solution at $t^* = 0.01 \ sec$ for probability to free-stream $p = 1$ (left) and $p = 0.5$ (right).

DownLoad: Full-Size Img PowerPoint

Figure 5.  Comparison of one realization of the MPC solution (RMPC without reactions/diffusion-only case) to the PDE solution at $t^* = 0.01 \ sec$ for probability to free-stream $p = 0.2$ (left) and $p = 0.1$ (right).

DownLoad: Full-Size Img PowerPoint

Figure 6.  Comparison of one realization of the MPC solution (RMPC without reactions/diffusion-only case) to the PDE solution at $t^* = 0.01 \ sec$ for probability to free-stream $p = 0.01$ for two different baths.

DownLoad: Full-Size Img PowerPoint

3.2.   Reaction-diffusion simulations

3.2.1.   Spatially inhomogeneous initial conditions

The number of cells in the $x$, $y$ and $z$ directions are $N_x = 200$, $N_y = 1$ and $N_z = 1$. The dimensionful time step is $\tau^* = 10^{-4} \ sec$, and the grid size is $h^* = 5\; \mu m$ with cubic cells. Fully reflective boundary conditions are imposed at all boundaries of the simulation volume and grid-shifting is implemented.

Using the best bath as reported in Table 1 for the diffusion-only case, the Schnakenberg reactions are added to the system dynamics involving two different species with different diffusion coefficients. The macroscopic reaction rates $k_1^*$ to $k_4^*$ in (3.2) are given the values $k_1^* = 4/\mu m^3 sec$, $k_2^* = 16/\mu m^3 sec$, $k_3^* = 6/sec$, and $k_4^* = 0.75 \mu m^6/sec$. This corresponds to rates for the number of particles in (3.4) as $k_1 = 500/sec$, $k_2 = 2000/sec$, $k_3 = 6/sec$, and $k_4 = 4.8 \hbox{x} 10^{-5} /sec$.

Figures 7 and 8 show a comparison of RMPC results with and without bath to the PDE solution with bath density chosen according to Table 1. For the first of the two figures, $D_U = 5 \ mm^2/sec$ and $D_V = 1 \ mm^2/sec$. In the second of the two figures, $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2/sec$. It is clear from the figures that the simulation results with bath lead to closer agreement with the PDE solution. However, the agreement is not nearly as good as in the non-reactive case. This could be a result of different local reactive dynamics, that are perhaps simply offset in time.

Figure 7.  Comparison of the RMPC solution (reaction-diffusion case) with and without bath to the PDE solution at $t^* = 0.01 \ sec$ for $U$ (left) with $D_U = 5 \ mm^2/sec$, and $V$ (right) with $D_V = 1 \ mm^2/sec$. Bath density is chosen as per Table 1.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Comparison of the RMPC solution (reaction-diffusion case) with and without bath to the PDE solution at $t^* = 0.01 \ sec$ for $U$ (left) with $D_U = 1 \ mm^2/sec$, and $V$ (right) with $D_V = 0.1 \ mm^2/sec$. Bath density is chosen as per Table 1.

DownLoad: Full-Size Img PowerPoint

To explore the differences in local reactive dynamics between the PDE and the stochastic results, Figures 9 and 10 show the phase portraits at two different locations for the RMPC simulations, compared to the PDE system solution, as well as the non-diffusive ODE system subject to the same initial condition as at the respective locations. It is clear from the figures that there are diffusive effects on the ODE dynamics, and that the phase portraits are comparable.

Figure 9.  Phase portraits at $x^* = 0.5 \ mm$ (center of domain) with $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2/sec$ of the RMPC solution (left) with best bath according to Table 1, and the PDE solution and non-diffusive ODE solution (right) for the same diffusivities and reaction rates.

DownLoad: Full-Size Img PowerPoint

Figure 10.  Phase portraits at $x^* = h^*/2 = 0.0025 \ mm$ (location of center of first RMPC cell) with $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2/sec$ of the RMPC solution (left) with best bath according to Table 1, and the PDE solution and non-diffusive ODE solution (right) for the same diffusivities and reaction rates.

DownLoad: Full-Size Img PowerPoint

In contrast, when no bath particles are used, the resulting phase portraits are shown in Figure 11 for the RMPC dynamics.

Figure 11.  Phase portraits at $x^* = 0.0025 \ mm$ (left) and $x^* = 0.5 \ mm$ (right) with $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2/sec$ of the RMPC solution with no bath particles.

DownLoad: Full-Size Img PowerPoint

It is clear that there is noise in the simulations, but that the RMPC reactive dynamics is comparable to the PDE solution at the two different locations in the simulation domain.

Finally, Figures 12 and 13 show the $x$-direction velocity histogram compared to the equilibrium velocity distribution for each of the respective species at the intermediate time $t^* = 0.01 \ sec$ (Figure 12), and at a much later time $t^* = 10 \ sec$ (Figure 13) for $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2/sec$ for the RMPC dynamics. One can see from Figure 12 that the equilibrium velocity distribution is well maintained for both $U$ and $V$ at the intermediate time for the results presented in Figure 8. Furthermore, the momentum transfer from $V$ to $U$ through reactions takes place only at much later times when the spatially homogeneous steady-state has been reached, and – as such – had a minimal effect on the diffusive behaviour presented in the Figure 8 results.

Figure 12.  Comparison of velocity histogram from RMPC simulations compared to the equilibrium velocity distribution for $U$ (left) and $V$ (right) at $t^* = 0.01 \ sec$ with $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2\sec$ with best bath.

DownLoad: Full-Size Img PowerPoint

Figure 13.  Comparison of velocity histogram from RMPC simulations compared to the equilibrium velocity distribution for $U$ (left) and $V$ (right) at $t^* = 10 \ sec$ with $D_U = 1 \ mm^2/sec$ and $D_V = 0.1 \ mm^2\sec$ when the best bath of Table 1 is used.

DownLoad: Full-Size Img PowerPoint

3.2.2.   Turing patterns

For the Turing pattern simulations, we perform dimensionless simulations with the number of cells in the $x$, $y$ and $z$ directions being $N_x = 200$, $N_y = 1$ and $N_z = 1$ respectively. Introducing dimensionless quantities for time ($t^* = T t$ for some time scale $T$), and dimensionless location ($x^* = L x$ for some length scale $L$) leads to a dimensionless system of equations. The dimensionless time step is taken as $\tau = 1$, and the dimensionless grid as $h = 1$ for the length of the sides of the cubic cells. Fully reflective boundary conditions are imposed at all boundaries of the simulation volume and grid-shifting is implemented. In order for the reaction rates to correspond to probabilities, we took $k_1 = 0.05$, $k_2 = 0.2$ and $k_3 = 0.00075$. Species $U$ and $V$ free-stream with different probabilities, where $k_B T/m = 1$ for both species, The initial condition consisted of $350$ U particles in every cell, while there were $200$ particles of V in the left half of the domain, and $100$ in the right half. No bath was used in these simulations.

With reaction rate $k_4 = 4.8 \hbox{x} 10^{-9}$, the simulation results for one realization in the case when $U$ diffuses with probability $0.001$ and $V$ with probability $0.1$ lead to a pattern with four peaks/dips for $U$ and $V$ respectively (see Figure 14). In this case, the PDE solution likewise predicts a pattern, although with more peaks/dips than the RMPC solution produced. In the Figure, we compare the RMPC solution with the PDE for which $pD_u = 0.0005$ and $pD_v = 0.3$ (rather than $0.05$).

Figure 14.  Particle number distributions of $U$ (left) and $V$ (right) for a Turing pattern compared to the PDE solution with $pD_u = 0.0005$, $pD_v = 0.3$ and $k_4 = 4.8 \hbox{x} 10^{-9}$.

DownLoad: Full-Size Img PowerPoint

Changing the reaction rate $k_4$ to $4.8 \hbox{x} 10^{-8}$, the PDE solution predicts a spatially homogeneous steady-state, while the RMPC simulations lead to a Turing pattern. Figure 15 shows the RMPC simulation compared to the PDE solution.

Figure 15.  Particle number distributions of $U$ (left) and $V$ (right) for a Turing pattern compared to the PDE solution with $pD_u = 0.0005$, $pD_v = 0.05$ and $k_4 = 4.8 \hbox{x} 10^{-8}$.

DownLoad: Full-Size Img PowerPoint

4.   Summary and Conclusions

Stochastic simulations using the reactive multiparticle collision (RMPC) dynamics were performed for a reaction-diffusion system obeying the Schnakenberg reaction dynamics, and subjected to a spatially inhomogeneous initial condition. An efficient global number density simulation approach was proposed so as to simulate a constant diffusion coefficient throughout the simulation domain. Constant diffusivities may be required for some applications where the diffusivity of the chemical species is known a priori, and additionally serves as a benchmark problem to compare with. The numerical approach was found to lead to correct diffusive motion of the particles provided that the mean free bath was large enough. For smaller mean free paths, correct diffusive motion was shown to be possible through the addition of bath particles. Adding reactions to the dynamics was shown to lead to correct phase portraits in key parts of the domain, although there was noise and time delays in approaching the homogeneous steady state of the system. Benefits of the global number density approach proposed here include numerical efficiency: there is one collision angle with random choice between plus/minus that angle locally, but the same choice of angles for every collision cell in the volume rather than a local method for which a collision angle is computed for each cell independently. Additionally, a novel means to simulate a chemical system for which the diffusivity of one species is significantly larger than the other was introduced through a probability to free-stream approach. The advantage of this approach, rather than varying $k_B T/m$, is that the same $k_B T/m$ value can be used for both chemical species in the system, and reactions do not change the Maxwell velocity distribution for any of the species in the system over long simulation times. This approach was used to successfully generate Turing patterns both when the PDE predicts a pattern, and when the PDE predicts a constant, spatially homogeneous solution.

The significance of the results presented here are as follows. Typically, stochastic simulation results are used so as to explore the breakdown of the PDE or ODE assumption. Typically any deviation is interpreted as noise-induced effects. Spatially extended models, or more specifically ones with position and velocity-dependence in the particle dynamics are scarce, with RMPC being a new and evolving simulation model for reaction-diffusion systems. The diffusive nature of RMPC in spatially inhomogeneous cases, which can naturally arise through reactive events, is a relatively unexplored area. The results presented here serve as a base for distinguishing between noise-induced differences in RMPC simulation results compared to PDE/ODE solutions, and spatially-dependent diffusion coefficients that RMPC dynamics might naturally produce. Additionally, this global approach has the advantage of avoiding any singularities that might arise in a local approach, in addition to being more computationally advantageous. Furthermore, for the first time, RMPC has been used to produce Turing patterns in a long-term simulation. Thus, RMPC can serve as a feasible simulation tool for both short and long-term simulations, and for a fairly large range of diffusivities that can be kept constant in a reaction-diffusion scenario. An advantage of this method includes more detailed system information compared to other approaches, with the ability to assess slip versus no-slip effects on the dynamics in a finite volume where boundary effects may have an impact in the overall dynamics, and the ability to model effects in cases where the mean-field PDE dynamics are no longer applicable.

Acknowledgments

The authors wish to thank the Natural Sciences and Engineering Research Council of Canada (NSERC; RGPIN-2017-04672), the Ryerson University Mathematics Department, and the Ryerson University Faculty of Science for funding support. Additionally, the simulations were made possible through the RAMLab (Ryerson Applied Mathematics Lab) thanks to an NSERC-RTI grant, as well as Faculty of Science funds for recent upgrades.

Conflict of interest

The authors declare that there is no conflict of interest in this paper.