An integration factor method for stochastic and stiff reaction–diffusion systems

Stochastic effects are often present in the biochemical systems involving reactions and diffusions. When the reactions are stiff, existing numerical methods for stochastic reaction diffusion equations require either very small time steps for any explicit schemes or solving large nonlinear systems at each time step for the implicit schemes. Here we present a class of semi-implicit integration factor methods that treat the diffusion term exactly and reaction implicitly for a system of stochastic reaction-diffusion equations. Our linear stability analysis shows the advantage of such methods for both small and large amplitudes of noise. Direct use of the method to solving several linear and nonlinear stochastic reaction-diffusion equations demonstrates good accuracy, efficiency, and stability properties. This new class of methods, which are easy to implement, will have broader applications in solving stochastic reaction-diffusion equations arising from models in biology and physical sciences.


Introduction
Complex patterns can be extensively found in nature, from the skin of zebrafish to the disposition of feather buds in chicks and hair follicles in mice.Often, those patterns are created by biochemical reactions along with diffusions of the molecules in a cellular or multi-cellular systems [1].Such biological systems, which may be described in reactiondiffusion equations, are constantly subjected to stochastic effects such as noises and environmental perturbations.The stochastic effects on the biochemical reactions at the single-cell level can result in heterogeneous responses of cellular populations and influence their behaviors [2].Previous studies on stochasticity reveal the adaptation of biological systems to noise, which can be characterized by the systems' strategies to combat noise, whether by attenuating or exploiting it [2] [3].For example, spatial stochastic effects help to either prompt the tight localization of proteins or enhance the response to the directional change of a moving pheromone input, resulting in a more robust cell polarization [4]; and the boundary of gene expression domains is sharpened as a result of gene-switching prompted by intracellular noise [5].It has become increasingly important to incorporate these stochastic effects into the reaction-diffusion equations for better understanding of biological systems.
One can describe a biological system in terms of the following stochastic reaction-diffusion equations: (1) where Ẇ (x, t) is a standard two-dimensional Wiener process.
One typical way of solving Eq.( 1) numerically is to apply central difference first to the diffusion and then use the temporal explicit scheme to solve the subsequent system [6], such as using the explicit Euler method [7] or the two, three, and four-stage explicit Runge Kutta schemes for the system containing additive noise and one-dimensional Wiener process in [8].Another common approach is using the Galerkin projection of the stochastic partial differential equation (SPDE) and then applying the numerical scheme to a finite-dimensional version of the SPDE.For example, Exponential Time Differencing (ETD) scheme may be applied to the Galerkin projection of the SPDE [9] [10].One other such example is the Lord and Rougemont scheme, which is derived through Galerkin projection and an integrating factor approach [11].This scheme is most effective for SPDE with Gevrey regularity, and more improvement may be made on such scheme by taking advantage of the Itô-Taylor expansion [11].
While explicit temporal schemes may be directly implemented for various spatial discretization, including finite element and Galerkin methods [12], to deal with the stability constraint associated with the diffusions, one can treat the diffusion term implicitly, while treating other terms explicitly [6] such as implicit Euler and Crank-Nicolson schemes [7].Higher order methods [13] can be achieved using Galerkin projection and the linear-implicit versions of strong Taylor schemes [10].Non-uniform time discretization on Brownian motion can also be obtained for implicit Euler scheme [14].
Stochastic stiffness arises from large differences in the magnitudes of Lyapunov exponents [15], resulting in the presence of different time scales.As in the deterministic case, explicit methods face step-size constraint when used to solve stiff stochastic systems [16].The timestep constraint can be improved with the modification of the stochastic term by adding more terms from the Itô-Taylor expansion for higher order of accuracy and stability.One of the most well-known schemes stemmed from this construction is the Milstein scheme [15].Treating the stochastic term implicitly is also one of the popular approaches [17] [18] [19] [20], albeit computationally expensive.Hence, a class of explicit methods known as Chebyshev methods are derived, which have better mean-square stability than explicit Euler method and are not as computationally expensive as implicit methods [21].A combination of different numerical schemes into one method can also be seen, such as the case of the Composite Euler scheme [22].For this scheme, at each temporal step, the stochastic differential equation is either solved by implicit Euler method or semi-implicit Euler method.The Composite method has similar order of convergence 1/2 to the Euler Maruyama method but better stability.
Most of the methods mentioned so far are derived to combat the stochastic stiffness through the improvement of the stochastic term, which can be costly and not as effective if the stiffness only occurs in the deterministic term.In such case, methods that treat the deterministic term implicitly while keeping an explicit treatment of the stochastic term are preferred [18].Here, we propose a new approach to the problem of stiffness caused by the deterministic term, more specifically the reaction term in Eq.( 1).The approach is based on the semi-Implicit Integration Factor (IIF) methods [23] [24] [25] [26], which has been found to be effective for the stiff reaction-diffusion equations with better stability constraints imposed on the time steps associated with both reaction and diffusion.In this approach, the time-step constraint for the diffusion term arising due to the inverse of the eigenvalues of the diffusion matrix, which can be large in magnitude, is resolved by treating the linear diffusion term exactly using Integration Factor (IF) methods.Such treatment results in an exponential function of the diffusion term and an integral of the nonlinear reaction term, which is then treated using implicit approach through the Lagrange interpolation to deal with its stiffness.Appropriate choices of approximation schemes lead to decoupling on the treatment for the diffusions and reactions such that one only needs to solve nonlinear systems with the size of the original PDEs.The IIF methods also have exceptional stability properties and its second-order version is absolutely stable.For higher-dimensional problems, the compact IIF (cIIF) method [24] is a great improvement on computational efficiency without altering the stability properties of the IIF methods [23].
In this paper, we exploit the simple structure of the IIF methods as well as their desirable stability properties and efficiency for solving the system in Eq. (1).Because of the nice decoupling properties in the IIF method, we will treat the deterministic diffusion and reaction terms in a similar fashion [23], while dealing with the stochastic term explicitly as in the Euler Maruyama method [27].We compare this approach with similarly constructed schemes whose main difference is in the treatment of the deterministic part of the equation, which can be approximated using ETD, Crank-Nicolson, or Implicit Euler methods.When all of the properties such as order of accuracy, mean errors, and stability region are taken into consideration, the new approach shows many advantages.We also take advantage of the low computational cost of the cIIF methods to similarly construct a stochastic method that can be applied to higher-dimensional problems.The paper is organized as followed.We first present the construction of the method for systems with additive noise or multiplicative noise, along with linear stability analysis and their comparisons with several other methods.Next, we compare the new method with other methods on linear SODEs and SPDEs on their orders of accuracy and stability constraints.Then, we use this approach to study a nonlinear activator-substrate system of two diffusion species and lastly, make our conclusion.

Construction of General Method
We consider the stiff reaction-diffusion equation with spatial white noise below: (2) where a∂ 2 U/∂x 2 is the diffusion term and a is a nonnegative constant, f(U) is the reaction term, and g(U)(∂ 2 W/∂x∂t) is the noise term of two possible forms: g(U)(∂ 2 W/∂x∂t) = σ(∂ 2 W/ ∂x∂t) for additive noise or g(U)(∂ 2 W/∂x∂t) = σU(∂ 2 W/∂x∂t) for multiplicative noise.Here, σ is a constant to describe the level of noise.Also, ∂ 2 W/∂x∂t denotes the mixed second-order derivative of the Brownian sheet.A one-dimensional Brownian sheet is a 2-parameter, centered Gaussian process B = B(s, t); s, t > 0 whose covariance is given by: (3) Before discussing the derivation of the numerical methods to solve Eq.( 2), we briefly review the Implicit Integration Factor methods (IIF) discussed in [23], which is crucial to the construction of the IIF methods for a stochastic system.Using the semi-discretized form dU/dt = aU + f(U) that is obtained after the discretization of the diffusion operator in space, we multiply both sides of the equation by the integrating factor e −at and integrating the equation over one time step from t n to t n+1 = t n + Δt to get (4) Using appropriate approximation of the integrands in one derives r thorder IIF scheme [23]: (5) with α n+1 ,α n , α n−1 ,…,α n−r+2 defined as (6) Similarly for the stochastic reaction-diffusion systems, we first discretize the space using m points with the spatial interval Δx.Let U t be a vector whose i th -entry is the value of the solution to Eq.(2) at the i th spatial point.A second-order central difference approximation of ∂ 2 U/∂x 2 in Eq.( 2) with periodic boundary condition U(x 0 , t) = U(x f , t) on the SPDE as in [28], where x 0 and x f indicate the endpoints of the spatial interval, leads to (7) where (8) Let and multiply both sides of this Eq.(7) by the integrating factor e −aMs .We then have (9) (10) (11) (12) Taking the integral of the left side gives (13) Letting Δt = t n+1 − t n and with some more simplification, the equation above becomes (14) All we have left is evaluating the right side of Eq.( 14).Observe the following numerical approximation of the noise part of Eq.( 14) [27] (15) To approximate the deterministic part of Eq.( 14), i.e , we apply the IIF strategy using Eqs.( 5) and (6).Coupling this evaluation with Eq.( 15), Eq. ( 14) becomes (16) where ΔW n = W t n+1 −W t n , U t n = U n , and α n+1 ,α n , α n−1 ,…,α n−r+2 are as described in Eq.( 6).
Let us denote ξ ̄n to be a standard normally-distributed random vector and n to be the indices of the temporal discretization points.We apply the standard Maruyama method to the noise term along with the first order IIF method, denoted as IIF1, or the second order IIF method, denoted as IIF2, to obtain IIF1-Maruyama method (17) IIF2-Maruyama method (18) When the stochastic integral in Eq.( 14) is approximated explicitly as in Eq.( 15), the strong order of convergence of the overall scheme is dominated by the root mean-square order of the increments ΔW n , which is one-half [15].For this reason, the strong order of convergence for both of our methods will be consistent with those of most other methods with the same approximation of 5 the stochastic term, i.e the Euler Maruyama method.Let us illustrate this through approximating an SODE of the similar form (19) using IIF1 and IIF2 methods respectively.By using the standard Maruyama approximation on the noise term, we obtain IIF1-Maruyama method: (20) IIF2-Maruyama method: (21) When the noise is additive, i.e. (22) the strong order of convergence for both methods is one, which is consistent with the strong order of convergence of the Euler Maruyama method [15].
When the noise is multiplicative, i.e. (23) both methods share the same order of convergence of one-half with the Euler Maruyama method [15].

2.2.1.
Additive noise-We first analyze the numerical stability of the IIF1-Maruyama and IIF2-Maruyama schemes in Eqs.( 20) and ( 21) when the stochastic differential equation has additive noise Eq. (22).
Following a previous study [15], we define to be a time discrete approximation of the solution U(t) with maximum step size δ > 0 starting at time t 0 at and to be the corresponding approximation starting at .Then is asymptotically numerically stable for a given stochastic differential equation if for any finite interval [t 0 , T] there exists a positive constant Δ a such that for each ε > 0 and δ ∈ (0, Δ a ) [15]: (24) and (25) with P(A) indicating the probability that event A occurs.We can analyze the asymptotic stability of a numerical stochastic scheme as we do for the A-stability of deterministic differential equations by studying the stability of the following class of complex-valued linear test equations [15]: (26) where λ is a complex number with ℝ(λ) < 0 and W is a real-valued standard Wiener process.
Suppose that a numerical scheme with equidistant step size Δt ≡ δ applied to test equation (26) with ℝ(λ) < 0 can be written recursively as: (27) where G is a mapping of complex plane ℂ into itself and are random variables that do not depend on for n = 0, 1, 2, …, then the set of complex values λΔt satisfying (28) is the region of absolute stability of that scheme [15].

Multiplicative
Noise-When the noise is multiplicative Eq.( 23), we analyze the stability of each method using meansquare stability analysis [15].A method is mean-square stable if lim n→∞ |U n | 2 ) = 0 where ․) denotes the expected value.To apply this technique to evaluating the stability region of both IIF-Maruyama methods aforementioned, we note that we can rewrite each method in the form: (30) Squaring and then taking expectation of both sides of Eq.( 30) coupled with the fact that W t is a standard Wiener process whose increment W(t) −W(s) is normally-distributed with mean 0 and variance t − s, we obtain (31) where Eq.( 31) demonstrates that lim n→∞ |U n | 2 ) = 0, i.e, the numerical method is mean-square stable if and only if H(a, b, σ, Δt) < 1 [17].
For the IIF1-Maruyama method, the mean-square stability condition becomes (32) For the IIF2-Maruyama method, the mean-square stability condition becomes (33) We plot the stability regions of both IIF-Maruyama schemes on a plane whose axes are aΔt and bΔt in Figure 1 and vary the value of σ 2 Δt.Note that the stability region in Figure 1 for each method is the region under the respective colored curve.The desired absolute stability region is the region where the diffusion and reaction coefficients are negative.In terms of Eq.( 19), this region is described as {(a, b) : a > 0 and b < 0}.In Figure 1A when there is no noise term, both methods are unconditionally stable with respect to this absolute stability region which is the inside of the square with dashed border.From Figure 1B and C, we observe that as the value of σ 2 Δt increases, the stability region of the IIF2-Maruyama method shrinks at a faster rate than the stability region of the IIF1-Maruyama method, resulting in the IIF1-Maruyama method having a larger stability region when the noise term is large enough.As a result, the IIF1-Maruyama method has a more desirable stability than the IIF2-Maruyama method in the case of more dominant noise.

Comparison with other methods in the case of Multiplicative Noise-For
the purpose of stability-region comparison, we present three other methods used to solve Eq. ( 19) and their constructions: The Euler Maruyama method [15] when it is applied to Eq.( 19) takes the form The order of accuracy for this method is 1/2 [15] and the mean-square stability analysis when noise is multiplicative, i.e Eq.( 23) gives the following stability condition: The next method is designed in a similar fashion to the construction of the IIF-Maruyama methods with a modification on the approximation of the deterministic integral term of Eq. ( 14).Direct application of the exponential time differencing method of order 2 on this term leads to Since the stochastic integral term is approximated explicitly as in the IIF-Maruyama methods, the order of accuracy for the overall scheme is 1/2.Mean-square stability analysis gives the following stability condition for the above method: (37) The last scheme mentioned here is constructed similarly to the Euler Maruyama method with the exception of the deterministic term being approximated using second-order Runge Kutta method.

RK2-Maruyama method
(38) The construction of the RK2-Maruyama scheme also exploits the explicit approximation of the stochastic term as that of the Euler Maruyama scheme, so the order of accuracy for the RK2-Maruyama scheme is still 1/2.The scheme's mean-square stability condition is: To illustrate the performance of the IIF-Maruyama schemes in terms of stability analysis in comparison with the above methods, we plot all the stability regions of each method for different values of σ 2 Δt on a plane whose axes are aΔt and bΔt (Figure 2).In the same figure, the region where unconditional stability is achieved for an ideal method is the region inside the box with dashed boundary.Figure 3 is the enlarged version of Figure 2 so we can observe better the changes in the absolute stability region for each method at different values of σ 2 Δt.
In Figure 2A, i.e when there is no noise term, only the IIF-Maruyama methods are unconditionally stable, which is consistent with the stability of the deterministic IIF methods.For very positive values of the diffusion term and very negative values of the reaction terms (in terms of Eq.( 19), this means that both a < 0 and b < 0), the Euler Maruyama, RK2-Maruyama, and ETD2-Maruyama methods achieve better stability than the IIF-Maruyama methods, as seen in the bottom left corner of each subplot of Figure 2A.However, the overall size of the absolute stability regions of the IIF-Maruyama methods is still larger than those of the other methods.
With the increasing size of the multiplicative noise term, the stability region for each method starts to shrink.More specifically, the stability regions for the Euler Maruyama and RK2-Maruyama methods start to shrink in width along the line b = −a (i.e along the direction in which the diffusion and reaction terms are equal) and become thin strips in Figure 2B and C. Both of these regions disappear completely in the next plot, i.e Figure 2D when the noise amplitude is high enough.
Meanwhile, the bottom-left corner of the stability region of the ETD2-Maruyama method recedes significantly as the noise amplitude increases, resulting in a greater loss of the absolute stability region than those of the IIF-Maruyama methods.Size comparison of the absolute stability regions from both Figure 2B-D and Figure 3B-D indicates that the stability region of the ETD2-Maruyama shrinks more than those of the IIF-Maruyama methods as σ 2 Δt = 1 increases from 0 to 1.
From Figure 2B-D, we observe that the IIF-Maruyama methods have the greatest region of absolute stability for any positive values of σ 2 Δt.Also, from Figure 3B-D, we observe that the stability regions for both IIF-Maruyama methods shrink at a slower rate than those of the other methods.As a result, both methods have the best absolute stability region for large noise amplitude, as demonstrated by Figure 2D, when σ 2 Δt = 1.In addition, the absolute stability region of the IIF1-Maruyama method shrinks more slowly than the IIF2-Maruyama, as evidenced by Figure 3B-D where the dashed box indicates the ideal absolute stability region.For this reason, the IIF1-Maruyama method achieves the largest absolute stability region for large noise amplitude out of all the methods.
We conclude that at different values of the noise term, the IIF-Maruyama methods outperform the other methods in terms of the region of absolute stability.In particular, at σ 2 Δt = 1, both the IIF-Maruyama methods have a much greater region of absolute stability than the rest of the methods.

Numerical Simulations
First, we compare the two IIF-Maruyama methods with the other methods when they are applied to Eq.( 2) for both cases of additive noise and multiplicative noise.Through choosing different values of a, which corresponds to the size of diffusion, and different values of b, which corresponds to the strength and stiffness of reactions, we evaluate the convergence and stability of IIF-Maruyama methods.

Tests on Stochastic Ordinary Differential Equations
Here, we implement various methods to solve the linear stochastic ODE Eq.( 19).
Comparisons will be made between the two IIF methods, the Euler Maruyama, and the ETD2-Maruyama methods.The comparisons concern the accuracy of these methods in situations where the degree of stiffness is high or the noise amplitude is great.All the simulations are done over 1000 independent paths with a time frame from 0 to 1 unless specified otherwise.Numerical experiments were needed in order to decide on a sufficiently large number of Brownian paths that will yield the desirable orders of convergence.All the results obtained in this sub-section remain consistent for a greater number of Brownian paths.This notion is confirmed when we increase the number of paths from 1000 to 2000 and subsequently 10000.

Additive Noise-Denote U
Δt to be the solution obtained numerically from using time step Δt.The order of convergence for additive noise is the value γ such that there exists a constant C where (40) for Δt sufficiently small.The order of convergence for the case, in which the explicit solution (e.g to Eq.( 19) with additive noise Eq.( 22)) is unknown, is estimated by In our simulations, we start with Δt = 2 −5 and decrease Δt by half for a total of 6 times.1000 independent Brownian paths are generated and the final solution U Δt on each path for each time step Δt is calculated.
Next, we study the accuracy and stablity for both IIF-Maruyama methods and compare them with the Euler Maruyama and ETD2-Maruyama methods in different scenarios, especially, in the case in which the reaction term is dominant and the system becomes stiff.
In Figure 4, we plot all the mean errors of the numerical solutions obtained from the IIF-Maruyama, Euler Maruyama, and ETD2-Maruyama methods while using different time steps in the scenario where the reaction term is heavily stiff, i.e when the magnitude of the reaction term is relatively large compared to the magnitudes of the diffusion and noise terms.Here, the mean error is defined to be U Δt − U Δt/2 | where U Δt is the numerical solution resulted from using each of the above methods with the time step Δt.From this figure, we observe that both IIF-Maruyama methods maintain a low mean error as the time step Δt increases in size.When Δt becomes too large, i.e when Δt = 1/4, the mean errors of the solutions from using Euler Maruyama and ETD2-Maruyama methods explode out of reasonable bounds.Meanwhile, at the same time step, the mean errors of the numerical solutions resulted from the two IIF-Maruyama methods remain consistently small when larger step size Δt is used.This figure demonstrates that the IIF-Maruyama methods are highly effective whenever the reaction-diffusion system has a dominant reaction term.
Figure 5 shows the orders of convergence for IIF1-Maruyama, IIF2-Maruyama, Euler Maruyama, and ETD2-Maruyama methods in different scenarios.Figure 5A and B represent the scenario where the noise amplitude is great compared to the magnitudes of the diffusion and reaction terms.Figure 5C and D represent the scenario where the magnitude of diffusion term is relatively large compared to those of the reaction and noise terms.In both scenarios, all the methods share an order of convergence of one as expected and no single method outperform the others.Figure 5 shows that both the IIF-Maruyama methods are comparable to other methods in terms of accuracy in the additive-noise case.

Multiplicative
Noise-For SODE with multiplicative noise, because the explicit solution for the linear stochastic ODE is known, the strong order of convergence can be estimated by the value γ if there exists a constant C such that: (42) for any fixed τ = nΔt ∈ [0, T] where T is the final time and for Δt sufficiently small.Let T = LΔt for some time step Δt [27].At τ = T, the order of convergence is calculated as the following First, we test the accuracy of the IIF-Maruyama, Euler Maruyama, and ETD2-Maruyama methods when the magnitude of the reaction term is large and plot the mean errors of the numerical solutions for each method in Figure 6A.The mean error for the multiplicative case is defined to be U n − U(τ)| from Eq.( 42).As in the additive case, when the reaction term is highly dominant, the IIF-Maruyama methods maintain low mean errors even when the time step Δt is relatively large.From this figure, we observe that when Δt is large enough, i.e Δt = 1/8, the mean errors of the solutions obtained from using Euler Maruyama and ETD2-Maruyama methods assume unreasonably large values and these two methods become unstable.At the same time step, both the IIF-Maruyama methods still maintain stability, as evidenced by the reasonable mean errors resulted from the numerical solutions.
We conclude that the IIF-Maruyama methods give reliable results in the case where the reaction-diffusion system is highly stiff in the reaction term.
Figure 6B shows that the orders of convergence for the IIF-Maruyama, Euler Maruyama, and ETD2-Maruyama methods are consistently one-half when the noise term is dominant.
There is no real advantage of choosing one method over another in this scenario.
When the magnitude of the diffusion term is relatively large compared to those of the reaction and the noise terms, we notice that the mean errors U n −U(τ)| obtained from the IIF-Maruyama methods take much smaller values than those of the Euler Maruyama and ETD2-Maruyama methods, as evidenced by Figure 6C and D. Although all the methods have relatively small mean errors, the IIF-Maruyama methods have the smallest mean error values and therefore are more accurate than the other methods.

Tests on Stochastic Partial Differential Equations
Now we apply IIF-Maruyama methods to the following stochastic PDEs and compare the IIF methods with two other methods. (44) where 0 ≤ x ≤ 1 and t ∈ [0, 0.125] along with a periodic boundary condition U(0, t) = U(1, t).
We compare the orders of convergence from solving Eq.( 44) in different scenarios among the following schemes: First-Order IIF-Maruyama method ( 17), Second-Order IIF-Maruyama method (18), Implicit-Euler Maruyama method, and Crank-Nicolson Maruyama method.The Implicit-Euler Maruyama method is constructed in the same way as the Euler Maruyama method with the deterministic term being approximated using the backward-Euler method of order one.The overall order of convergence for this method is 1/2 due to the explicit approximation of the stochastic term as in the case of the Euler Maruyama method.Letting M denote the diffusion matrix after applying finite difference to Eq.( 44), the method takes the following form, The construction of Crank-Nicolson method is also similar to that of the Euler Maruyama method, with the exception of the approximation of the deterministic term using the Crank-Nicolson method.For the same reason as the implicit-Euler Maruyama method, the overall order of convergence remains 1/2 for this method: We do not show the numerical results of the explicit Euler Maruyama scheme Eq.( 34) and ETD2-Maruyama scheme Eq.(36) due to their disadvantages in stability and the associated computational cost.
To compute the order of convergence for each scheme mentioned, we use five different values of the number of spatial steps: N 1 , N 2 ,…N 5 where N 1 is a power of 2 and N i+1 = 2N i for i = 1, …, 4, and let the time step Δt = 1/(4N i ).The solutions are numerically calculated over the time frame [t 0 , t f ] and generated over m different realizations of the Brownian sheet W(t, x).More information on how the Brownian sheet is generated can be found in [28].We record the difference at the final time between two solutions obtained using N i and N i+1 spatial steps and store this difference under the variable S i for i = 1, …4.The difference is the sum over m realizations of the sum of squared differences of the approximated solutions over N 1 spatial points, which are common to all solutions.Therefore, we obtain (47) which offers a mean to calculate the numerical error of the scheme.Note that indicates the approximated solution at space step x k = k/N 1 and at final time, which is obtained from using N i spatial steps and j th independent realization of the Brownian sheet [7].Then the order of convergence can be estimated by log 2 R where the ratio R = S i /S i+1 .When using this method of computing the order of convergence, both Implicit-Euler Maruyama and Crank-Nicolson Maruyama schemes converge with an order of 1/2 for both additive and multiplicative noises [7].As a result, we will use 1/2 as the standard value of the order of convergence in the subsequent numerical comparisons.In this sub-section, we fix m = 100.From our experimentation, increasing the value of m has no effects on the orders of the convergence of each scheme.However, the values of will increase since these quantities depend on the value of m.In our tests, when m = 500, the values of are roughly five times larger than those obtained with m = 100.Similarly, if we increase m to 1000, are about ten times larger than their corresponding values when m = 100.
To compare the orders of convergence, we observe the following scenarios with both noises: how the order of convergence for each scheme is affected when the degree of stiffness increases, and whether each method still performs satisfactorily with large noise amplitude.7A and B plots the values obtained from the numerical solutions of Eq.( 44) with both additive and multiplicative noises versus the size of the space step when the reaction term is stiff with respect to the diffusion and noise terms.The numerical methods that are applied here are the IIF1-Maruyama, IIF2-Maruyama, Implicit-Euler Maruyama, and Crank-Nicolson Maruyama methods.For both types of noise, all methods converge with a rate of approximately 1/2.The Implicit-Euler Maruyama scheme along with the two IIF-Maruyama schemes have an advantage over the Crank-Nicolson Maruyama scheme when the values are considered.

Stiff reaction-Figure
In the interest of demonstrating the effectiveness of the IIF-Maruyama methods, we compare the orders of convergence between all of the aforementioned methods Eqs.( 17), ( 18), (45), and (46) when the reaction term is extremely stiff for both additive and multiplicative noises.When noise is additive, all methods have an order of convergence of 1/2.However, the Implicit-Euler Maruyama and the two IIF-Maruyama methods perform much better than the Crank-Nicolson Maruyama method in terms of the mean errors , as seen in Figure 7C.In Figure 7D, for multiplicative noise, the Implicit-Euler Maruyama and the IIF1-Maruyama methods converge at a much faster rate than 1/2 while the Crank-Nicolson Maruyama and the IIF2-Maruyama schemes maintain the 1/2 order of convergence.When using fewer number of space steps, the values S i obtained from the Implicit-Euler Maruyama and IIF1-Maruyama methods are not as good as those of the other methods due to the large convergence rate.For example, when the spatial step sizes are 1/64, 1/128, and 1/256, the values S 1 and S 2 obtained from the IIF1-Maruyama and Implicit-Euler Maruyama methods are larger than those of the IIF2-Maruyama and Crank-Nicolson Maruyama methods.On the other hand, the IIF2-Maruyama method consistently has the smallest values for , making it the most desirable method for solving a stochastic partial differential equation with an extremely stiff reaction term whose noise term can be either additive or multiplicative.

Strong diffusion-
We obtain the numerical solutions for Eq.( 44) using the four methods when the diffusion coefficient is large.Then we compare the values from Eq.( 47) and the orders of convergence of all methods by plotting versus the size of the space step along with a reference line of slope one-half in Figure 8A and B. The Crank-Nicolson Maruyama method is slightly more stable than the other methods for some large values of spatial step size but does not maintain this stability if the space step assumes a larger value than those shown in this figure.Both the IIF-Maruyama methods and the Implicit-Euler Maruyama method achieve the best values for , while the Crank-Nicolson Maruyama method has significantly larger compared to them.When the spatial step size assumes a small-enough value, the IIF-Maruyama schemes and the Implicit-Euler Maruyama schemes have comparable orders of convergence with that of the Crank-Nicolson Maruyama scheme.With both the mean errors and orders of convergence taken into consideration, it is more advantageous to choose the IIF-Maruyama methods and the Implicit-Euler Maruyama method over the Crank-Nicolson Maruyama scheme.8C and D contains similar plots to Figure 8A and B in the case where Eq.( 44) has a large noise term.For both additive and multiplicative noises, all methods have an order of convergence of one-half.Also, for all methods, the values of are slightly larger than the corresponding values obtained when the reaction term or diffusion term is stiff.Since the calculation of contains a double sum, the magnitude of could become quite large.Taking this into consideration, when the noise term has large magnitude, the IIF-Maruyama methods and the Implicit-Euler Maruyama method outperform the Crank-Nicolson Maruyama method significantly in terms of mean errors and thus are preferred.

A Turing patterning system with noise 4.1. 1-dimensional activator-substrate system
Finally we apply the IIF2-Maruyama scheme to a Turing patterning system that contains noise to study how noise may affect the formation of patterns.We use the activator-substrate system as an example, whose non-dimensional form is as followed [29] [30]: (48) (49) state counterparts.Unlike the previous study for the deterministic equations [29] in which fluctuations in the initial conditions are critically important in generating the patterns, we fix the initial conditions A 0 = A * and S 0 = S * defined in Eq.( 50).
Interestingly, we obtain similar patterns in spite of the uniform initial conditions when the relatively small values to the noise coefficients, i.e ε A = 0.005 and ε S = 0.01, are given.The six different combinations of patterns for the long time solutions A and S that exist can be seen in Figure 9.We note that these six different combinations of patterns are the same inhomogeneous steady state patterns obtained from solving the deterministic equations Eqs. ( 48) and (49) [29].Hence, adding multiplicative noise to the activator-substrate system enables us to obtain the inhomogeneous steady-state patterns that are otherwise obtained through the fluctuations of the initial conditions, as previously predicted [31].
Because each stochastic solution may reach a different steady state even with the same initial conditions in a deterministic form, we perform 100, 500, and 1000 simulations to see which combination of patterns shows up more frequently.In Table 1, we fix the initial conditions as in Eq.( 50) and choose the spatial step size to be dx = 10/2 6 and dx = 10/2 7 .This change in dx does not affect the frequency of presence of each combination of patterns.For Table 2, the values of the initial conditions are randomly permuted as in Eq.( 52) for each independent path, while dx is fixed at 10/2 7 .From these two tables, we see that the first combination of patterns is consistently the most favored type of patterns, with a > 30% chance of occurrence, with the second combination of patterns being the second most typical combination of patterns.In addition, the frequency of appearance of each type of patterns is independent of the effects of extra fluctuations on the initial conditions, as evidenced by Table 2.This implies that the first combination of patterns is likely to make up the standard type of patterns that the long-term activator and substrate solutions are supposed to assume.The lack of robustness in pattern formation of the activator and substrate levels can be improved by adding growth factor to the system, in particular apical growth in the case of intrinsic noise [31].As the domain grows, the space between the activated regions (characterized by the activator maxima or the substrate minima) is enlarged while the substrate concentration is quickly diffusing and increasing.This increase in substrate prompts a higher production of the activator at the side of the maxima in comparison to its center, resulting in the movement of the activator maxima to regions with higher substrate concentration [30].Attributing the appropriate type of growth to the system can help stabilize the pattern formation over time where the robust patterns for both activator and substrate levels are the first combination of patterns.

4.1.2.
Computational efficiency-Next, we discuss the computational efficiency of the two IIF-Maruyama methods used to solve the Activator-Substrate system with multiplicative noise, Eqs.(54) and (55), by comparing their performances with that of an explicit method, which we choose to be the Euler Maruyama method.We keep all the parameter values as described previously Eq.(53) with the exception of time, which is changed to t ∈ [0, 1].In Table 3, the mean errors, the order of convergence, and the computational time for each method are recorded.To calculate the mean errors and the order of convergence, we carry out 100 different simulations and apply Eq.(47).We choose dt = dx/2 for the IIF1-Maruyama method, dt = dx for the IIF2-Maruyama method, and dt = (dx) 2 /3 for the Euler Maruyama method to ensure convergence.In this table, we denote N to be the number of spatial grid points that partition the interval (0,10).Besides N = 2 3 , all methods display similar orders of convergence and mean errors of similar magnitude.Next, we discuss the computational time in seconds of each method.Each time listed in our table is the total time each method takes to compute the solution over 100 different Brownian paths.For the computational time of both IIF-Maruyama methods, we include the time it takes to calculate the exponential matrix using the Matlab function expm.We observe that with coarser spatial grids, i.e when N = 2 3 and N = 2 4 , the Euler Maruyama method surpasses both IIF-Maruyama methods in terms of computational effort.With finer grids, i.e N = 2 6 and 2 7 , it takes the Euler Maruyama method twice as long as it takes the IIF2-Maruyama method in computing the solutions over 100 Brownian paths.Similarly, the IIF1-Maruyama method is put at a disadvantage with respect to computational speed when the grids are coarse but quickly catches up to the Euler Maruyama method with refined grids.We note that the second-order IIF-Maruyama method is more efficient than the first-order IIFMaruyama method.In addition, the second-order IIF-Maruyama method catches up to the Euler Maruyama method much faster in improving its computational speed, which is demonstrated by the similar speeds between these two methods when N = 2 4 .Meanwhile, the computational speed of the IIF1-Maruyama method does not catch up to that of the Euler Maruyama method until N = 2 5 .In short, due to the restriction of the temporal step size that is required to maintain numerical stability, the Euler Maruyama scheme is less computationally efficient than the IIF1-Maruyama and the IIF2-Maruyama schemes when a finer spatial grid is required.Between the two IIF-Maruyama methods, the second-order IIF-Maruyama method is more desirable for its computational efficiency than the first-order IIF-Maruyama method.

2-dimensional activator-substrate system
For two or three-dimensional systems, direct application of IIF is costly because the storage and computation of the exponential matrix e aΔt in the IIF methods may become very large.Similarly to solving the deterministic systems in two or three dimensions, here we use the compact integration factor methods (cIIF) [24], in which the discretized diffusion operator is represented in a compact form that requires storage only proportional to the number of unknowns instead of the square of the number of unknowns in the case of non-compact IIF methods for the exponentials of matrices.cIIF methods can be combined with the Maruyama method in the same manner as the integration of the IIF methods with the Maruyama method, and both the cIIF and the IIF methods share the same desirable stability properties [24].
Here we construct the cIIF-Maruyama methods by estimating the deterministic diffusion and reaction terms using the cIIF methods and the stochastic term using the explicit Maruyama approximation.To demonstrate the efficiency of the cIIF-Maruyama methods, we apply the cIIF2-Maruyama method to the two-dimensional version of the activator-substrate system with no-flux boundary conditions Eqs.(54) and (55) presented in the previous section.Similarly to the one-dimensional case, we compute the solutions over the space (0, 10) × (0, 10) using the steady states in Eq.(50) as the initial conditions.The time window for the simulation is set to be t ∈ [0, 200] and dt = dx where dx = dy = 2 −7 × 10.The rest of the parameter values in Eq.( 53) remain the same.Figure 10 displays the contour plot of one of the final patterns of the solutions A and S.
The computational time it takes for the cIIF2-Maruyama method to compute one stochastic solution to the two-dimensional activator-substrate system is 126.789seconds for this case.If we use the Euler-Maruyama method to solve this system using the same spatial step size, it takes 206.740 seconds due to the restriction on the temporal step size, which we set to be dt = (dx) 2 /5.In summary, due to the stability property, the cIIF-Maruyama methods are as efficient in solving two dimensional stochastic reaction-diffusion systems as the onedimensional systems.

Conclusions and Discussion
By convention, stochastic stiffness is defined to be the result of the different time scales caused by the large discrepancies in the magnitudes of the Lyapunov exponents [15] [18].For that reason, stochastic stiffness can occur in the deterministic term, stochastic term, or both.When solving a stochastic differential equation, the problem of stiffness that stems from the stochastic term has been studied previously [15] [16] [17] [18] [19] [20] [21] [22].
Here, we have focused on the treatment of stiffness of the reaction term for a stochastic reaction-diffusion system.By taking advantage of the existing semi-implicit integrating factor method that is both computationally efficient and absolutely stable at solving stiff deterministic reaction-diffusion systems, we have developed a new class of temporal schemes for reaction-diffusion systems with both additive and multiplicative noises.Similarly to the deterministic case, the new numerical schemes presented remove the restriction imposed on the temporal step size by the linear diffusion term by treating this term exactly while dealing with the stiff reaction term through an implicit approximation.Numerical comparisons show that the construction using the IIF technique to approximate the deterministic term allows the new methods to achieve better stability and good efficiency.While the explicit treatment of the diffusion in IIF naturally leads to good approximations on strong diffusion, the new IIF-Maruyama methods mainly offer an efficient approach to deal with stiff reactions in a reaction-diffusion systems.In general, this method is mostly effective when the reactions are very stiff while diffusion is still important in a stochastic reaction-diffusion system.
The approach used here in combining IIF for reaction and diffusion and Maruyama for the stochastic terms can be adapted in a straightforward fashion for compact IIF (cIIF) methods [24], which is effective for systems in two or three dimensions.With the compact representation for the differential operators, it would be more efficient in simulating 2D and 3D using cIIF than IIF-Maruyama.Another improvement on IIF-Maruyama is its order of convergence.
In the case of multiplicative noise, the order of convergence of the IIF-Maruyama methods is confined to a value of one-half.With the Euler Maruyama method, the order of convergence can be improved by adding more terms from the Ito-Taylor expansion to construct methods with higher order of accuracy such as the Milstein method [15].Integrating such approach with IIF method might lead to higher order of accuracy with similar stability property of IIF-Maruyama.However, because the diffusion term is not welldefined for the Brownian sheet due to a lack of a well-defined spatial derivative for the Wiener process, a direct application of the Milstein method to stochastic PDEs may only lead to half-order of convergence.To deal with this difficulty, one might need to use the Q-Wiener process instead of the Brownian sheet to approximate the diffusion [32].Similarly to the IIF methods, which can be used for the deterministic systems in various forms, spatial dimensions, and in combination with other approaches for treating additional terms such as convections, IIF-Maruyama methods may have broad applications in simulating stochastic partial differential equations in various forms and containing stiff reactions.The stability regions of both IIF-Maruyama methods described in Eqs.( 17) and ( 18) for multiplicative noise.The stability region lies below the corresponding colored curve.The desired absolute stability region is the region inside the square with dashed-border.Comparison of stability regions of the following methods: IIF1-Maruyama, IIF2-Maruyama, Euler Maruyama, RK2-Maruyama, and ETD2-Maruyama used to solve Eq.( 19) with multiplicative noise.The stability region for each method is shaded blue.The ideal absolute stability region is the region inside the dashed box.One of the final patterns obtained for long-term solutions A and S when solving the twodimensional stochastic system in Eqs.( 54) and ( 55).The values for the noise coefficients are ε S = 0.01 and ε A = 0.005.

Table 3
The mean errors S defined in Eq.(47), orders of convergence, and computational cost when solving the activator-substrate system with noise described in Eqs.(54) and (55) obtained by each method: IIF1-Maruyama, IIF2-Maruyama, and Euler Maruyama.The results are computed over 100 Brownian paths with N denoting the number of spatial grid points that partition the interval (0,10).

IIF2-Maruyama
denote U L,Δt to be the numerical solution at t = T using time steps of size Δt each and U L,Δt/2 the numerical solution at time T using time steps of size Δt/2 each.To acquire the value for U L −U(T)|, we take the mean of |U L −U(T)| over 1000 independent Brownian paths, hence we call U L − U(T)| the mean error of the numerical solution.

Figure 3 .
Figure 3.Enlarged version of Figure2.In this figure, the stability regions are enlarged for the following methods: IIF1-Maruyama, IIF2-Maruyama, and ETD2-Maruyama.The stability region for each method is shaded blue.The ideal absolute stability region is the region inside the dashed box.

Figure 4 .
Figure 4. Comparison of the mean errors U Δt −U Δt/2 | of the numerical solutions to Eq.(19) with additive noise obtained from the following methods: IIF1-Maruyama, IIF2-Maruyama, Euler Maruyama, and ETD2-Maruyama, when the reaction term is heavily stiff.The inserted figure shows the mean error comparison between different methods in more detail.The parameter values are as followed: a = 1, b = −10, and σ = 0.1.

Figure 5 .
Figure 5.Comparison of the orders of convergence of the following methods: IIF1-Maruyama, IIF2-Maruyama, Euler Maruyama, and ETD2-Maruyama used to solve Eq.(19) with additive noise.Subplots A and B represent the orders of convergence of all the methods in the scenario that the noise amplitude is relatively large compared to the magnitudes of the diffusion and reaction terms, whose values are fixed to be a = 0.1 and b = −1.Subplots C and D represent the scenario where the magnitude of the diffusion term is great compared to those of the reaction and noise terms, whose values are fixed to be σ = 0.1 and b = −1.

Figure 6 .
Figure 6.Comparison of the mean errors and orders of convergence between the following methods: IIF1-Maruyama, IIF2-Maruyama, Euler Maruyama, and ETD2-Maruyama, used to solve Eq. (19) with multiplicative noise.Subplot A is the comparison of the mean errors U n − U(τ)| of the numerical solutions when the reaction term is heavily stiff.The parameters used here are a = 1, b = −20, and σ = 0.1.In the scenario that the noise amplitude is relatively large compared to the magnitudes of the diffusion and reaction term, all methods display similar orders of convergence, as seen in subplot B. The parameter values for this subplot are a = 0.1, b = −0.02,and σ = 1.Finally, subplots C and D compare the mean errors U n − U(τ)| of the numerical solutions when the diffusion term is dominant using fixed parameters σ = 0.1 and b = −2.The time span for the simulations in subplot A is one and the time span used in subplots B-D is one-half.For subplots A and D, the inserted images show the mean errors of each method in more detail.

Figure 7 .
Figure 7.Comparison of the values from Eq.(47) and orders of convergence of the following methods: IIF1-Maruyama, IIF2-Maruyama, Implicit-Euler Maruyama, and Crank-Nicolson Maruyama used to solve Eq.(44) with both additive and multiplicative noises in the scenario where the reaction term is heavily stiff.Subplots A and B show the plots of of all the methods when the reaction term b = Subplots C and D display plots of of all the methods when b = −50.For all the subplots, the values of the reaction and noise terms are fixed to be a = 1 and σ = 0.1.In subplots A and B, the time frame is chosen to be t ∈ [0, 0.125] and for the remaining two subplots, t ∈ [0, 0.025].Also, each plot contains the reference line of slope one-half for the purpose of order of convergence comparison.

Figure 8 .
Figure 8.Comparison of the values from Eq.(47) and orders of convergence of the following methods: IIF1-Maruyama, IIF2-Maruyama, Implicit-Euler Maruyama, and Crank-Nicolson Maruyama used to solve Eq.(44).Subplots A and B are the plots of of all the methods when the diffusion term is stiff.The values of the diffusion, reaction, and noise terms are fixed to be a = 20, b = −1, and σ = 0.1 for all subplots.Subplots C and D display plots of of all the methods when the noise term assumes a large value.For these subplots, the values of the diffusion, reaction, and noise terms are fixed to be a = 2, b = −1, and σ = 1.In this figure, all the simulations are run for 0.125 time units.Also, each plot contains the reference line of slope one-half for the purpose of order of convergence comparison.

Figure 9 .
Figure 9.The six different combinations of steady-state patterns for long-term solutions A and S to the one-dimensional stochastic system in Eqs.(54) and (55).The values for the noise coefficients are ε S = 0.01 and ε A = 0.005.

J
Comput Phys.Author manuscript; available in PMC 2016 August 15.

Table 2
Percentage of occurrence of each combination patterns for 100, 500, and 1000 different simulations.In this table, dx = 2 −7 and the initial conditions are permuted as in Eq.(52).
J Comput Phys.Author manuscript; available in PMC 2016 August 15.