Stability-Optimized High Order Methods and Stiffness Detection for Pathwise Stiff
Stochastic Differential Equations
Chris Rackauckas

Qing Nie

Department of Mathematics
Massachusetts Institute of Technology
Cambridge, Massachusettss 02139
Email: crackauc@mit.edu

Department of Mathematics
University of California, Irvine
Irvine, California, 92697
Email: qnie@uci.edu

Abstract—Stochastic differential equations (SDE) often exhibit
large random transitions. This property, which we denote as
pathwise stiffness, causes transient bursts of stiffness which
limit the allowed step size for common fixed time step explicit
and drift-implicit integrators. Here we present a HPC-driven
method for deriving high strong order methods for stochastic
differential equations. Utilizing GPU-accelerated global optimization, we numerically solve a constrained optimization
problem which results stability-optimized adaptive methodsof strong order 1.5 for SDEs. The resulting explicit methods
are shown to exhibit substantially enlarged stability regions
which allows for them to solve pathwise stiff biological models
orders of magnitude more efficiently than previous methods
like SRIW1 and Euler-Maruyama. These methods are benchmarked on a range of semi-stiff problems and demonstrate
speedups between 6x previous adaptive algorithms while showing computationally infeasibility of fixed time step integrators
on some of these test equations.

1. Introduction
Stochastic differential equations (SDEs) are dynamic
equations of the form
dXt = f (t, Xt )dt + g(t, Xt )dWt ,where Xt is a d-dimensional vector, f : Rd → Rd is the drift
coefficient, and g : Rd → Rd×m is the diffusion coefficientwhich describes the amount and mixtures of the noise
process Wt which is a m-dimensional Brownian motion.
SDEs are of interest because they can exhibit behaviors
which are not found in deterministic models. For example,
an ODE model of a chemical reaction network may stay at
a constant steady state, but in the presence of randomness
the trajectories may be switching between various steadystates [1], [2], [3]. These types of models can capture crucial
features such as transitions to cancerous states [4], discoverrare cellular behaviors [5], and enhance personalized drug
metabolism predictions for precision healthcare [6].
In many cases, these unique features of stochastic models are pathwise-dependent and are thus not a property ofthe evolution of the mean trajectory. However, these same

effects cause random events of high numerical stiffness,
which we denote as pathwise stiffness, which can cause
difficulties for numerical integration methods.
A minimal example of pathwise stiffness is demonstrated
in the equation
dXt = [−1000Xt (1 − Xt ) (2 − Xt )] dt + g(t, Xt )dWt ,(1)
where X0 = 2, t ∈ [0, 5], and with additive noise g(t, Xt ) =10 where a sample trajectory is shown in Figure 1. Thisequation has two stable steady states, one at X = 0 andanother at X = 2, which the solution switches between
when the noise is sufficiently large. While near a steady state
the derivative is approximately zero making the problem
non-stiff, during these transitions the derivative of the drift
term reaches a maximum of ≈ 400. This means that inorder to be stable, explicit Stochastic Runge-Kutta (SRK)
must have a small ∆t. This display of large, transient,and random switching behavior in a given trajectory causes
stochastic bursts of numerical stiffness, a phenomena which
we will denote pathwise stiffness. The fixed time step
Euler-Maruyama method would require dt < 4 × 10−3 tobe stable for most trajectories, thus requiring greater than
2 × 104 steps to solve this 1-dimensional SDE. In manycases the switching behavior can be rare (due to smaller
amounts of noise) or can happen finitely many times likein the multiplicative noise version with g(t, Xt ) = 10Xt .
Yet even if these switches are only a small portion of
the total time, the stability requirement imposed by their
existence determines the possible stepsizes and thus has alarge contribution to the overall computational cost. Whileimplicit methods can be used to increase the stability range,
this can vastly increase the overall computational cost ofeach step, especially in the case large systems of SDEs likediscretizations of stochastic reaction-diffusion equations. In
addition, implicit solvers have in practice a smaller stability
region due to requiring convergence of the quasi-Newton
solvers for the implicit steps. This problem is mitigatedin ODE software by high-quality stage predictors given by
extrapolation algorithms for good initial conditions for the
Newton steps [7]. However, there are no known algorithms
for stage predictors in the presence of large noise bursts andPathwise Stiffness Example

that these adaptive stability-optimized SRK methods canadequately solve transiently stiff equations without losing
efficiency in non-stiff problems. Together we test on semistiff equations with 2 to 6 × 20 × 100 SDEs from biologicalliterature and show a speedup around 6x over the previous
adaptive SRIW1 algorithm, and demonstrate the infeasibilityof common explicit and implicit methods (Euler-Maruyama,
Runge-Kutta Milstein, Drift-Implicit Stochastic θ-Method,and Drift-Implicit θ Runge-Kutta Milstein) found as thebasis of many SDE solver packages [16], [17], [18].

2.0

Xt

1.5
1.0
0.5
0.0
0

1

2

t

3

4

5

Figure 1. Example of a Pathwise Stiff Solution. Depicted is a sample
trajectory of Equation 1 solved using the SOSRI methods developed inthis manuscript with reltol = abstol = 10−2 .

thus we will demonstrate that classic implicit solvers havea form of instability. Additionally, previous literature questioned the applicability of L-stable integrators to stochastic
differential equations due to high error in the slow variables[8]. Thus both fixed time step explicit and implicit solvers
are inadequate for efficiently handling this common class of
SDEs.
Since these features exist in the single trajectories of therandom processes, methods which attempt to account forthe presence of such bursts must do so on each individual
trajectory in order to be efficient. In previous work, theauthors have shown that by using adaptive time-stepping,a stochastic reaction network of 19 reactants is able tobe solved with an average time step 100,000 times larger
than the value that was found necessary for stability during
the random stiff events for a high order SRK method [9].
This demonstrated that the key to solving these equations
efficiently required controlling the time steps in a pathwisemanner. However, the methods were still largely stabilitybound, meaning the chosen tolerances to solve the model
were determined by what was necessary for stability andwas far below the error necessary for the application. The
purpose of this investigation is to develop numerical methods with the ability to better handle pathwise stiffness andallow for efficient solving of large Monte Carlo experiments
without resorting to drift-implicit methods.
We approach this problem by developing adaptive
stability-optimized SRK methods with enlarged stability
regions. This builds off of similar work for ODE integrators
which optimize the coefficients of a Butcher tableau to giveenhanced stability [10], [11], [12]. Similar to the RungeKutta Chebyschev methods [7] (and the S-ROCK extension
to the stochastic case [13], [14], [15]), these methods aredesigned to be efficient for equations which display stiffness
without fully committing to implicit solvers. Given the complexity of the stochastic stability equations and order conditions, we develop a novel and scalable mechanism for thederivation of “optimal” Runge-Kutta methods. We use thismethod to design stability-optimized methods for diagonal
noise SDEs. We show through computational experiments

2. Adaptive Strong Order 1.0/1.5 SRK Methods for Diagonal Noise SDEs
The class of methods we wish to study are the adaptive
strong order 1.5 SRK methods for diagonal noise [9], [19].
Diagonal noise is the case where the diffusion term g is diagonal matrix σi Xti and includes phenomenological noisemodels like multiplicative and affine noise. The diagonal
noise methods utilize the same general form and orderconditions as the methods for scalar noise so we use theirnotation for simplicity. The strong order 1.5 methods for
scalar noise are of the form
s


X
(0)
(0)
Xn+1 = Xn +αi f tn + ci h, Hi
+
(2)
i=1


s 
X
(4) I(1,1,1)(3) I(1,0)(1)(2) I(1,1)+ βiβi I(1) + βi √ + βi
h
h
h
i=1
(3)


(1)
× g t n + ci h(4)with stages(0)Hi= Xn +

s
X



(0)
(0)
(0)
Aij f tn + cj h, Hj
h+

(5)

j=1
s
X



(0)
(1)
(1) I(1,0)
Bij g tn + cj h, Hj
h
j=1

(1)

Hi

= Xn +

s
X



(1)
(0)
(0)
Aij f tn + cj h, Hj
h+

(6)
(7)

j=1
s
X


√
(1)
(1)
(1)
h
Bij g tn + cj h, Hj

j=1

where the Ij are the Wiktorsson approximations to theiterated stochastic integrals [20]. The tuple of coefficients
A(j) , B (j) , β (j) , α thus fully determines the SRK method.
These coefficients must satisfy the constraint equations described in [9], [19] in order to receive strong order 1.5.
These methods are appended with error estimates


X
(0)
(0)
ED = ∆t(−1)σ(i) f tn + ci ∆t, Hi
i∈I1

EN

=

 

X  (3) I(1,0)(4) I(1,1,1)
(1)
(1)
+ βi
g tn + ci ∆t, Hi
βi
∆t
∆t

i∈I2

and the rejection sampling with memory (RSwM) algorithm
to give it fully adaptive time-stepping [9]. Thus unlike inthe theory of ordinary differential equations [21], [22], [23],[24], [25], the choice of coefficients for SRK methods does
not require explicitly finding an embedded method whendeveloping an adaptive SRK method and we will therefore
take for granted that each of the derived methods is adaptive.

Notice that for N, M → ∞, f is the area of the stability
region. Thus we define the stability-optimized diagonal SRK
method as the set of coefficients which achieves

3. Optimized-Stability Order 1.5 SRK Methods with Diagonal Noise

However, we impose a few extra constraints to add robustness to the error estimator. In all cases we impose(0) (1)(0)0 < ci , ci < 1 . Additionally we can prescribe c4 =(1)
c4 = 1 which we call the End-C Constraint. Lastly, we can
(j)
(j)
(j)
(j)
prescribe the ordering constraint c1 < c2 < c3 < c4
which we denote as the Inequality-C Constraint.
The resulting problem is a nonlinear programming problem with 44 variables and 42-48 constraint equations. The
objective function is the two-dimensional integral of a discontinuous function which is determined by a polynomialof in z and w with approximately 3 million coefficients.
To numerically approximate this function, we calculated
the characteristic function on a grid with even spacing dx
using a CUDA kernel and found numerical solutions to the
optimization problem using the JuMP framework [27] with
the NLopt backend [28]. A mixed approach using manysolutions of the semi-local optimizer LN AUGLAG EQ[29], [30] and fewer solutions from the global optimizer
GN ISRES [31] were used to approximate the optimality ofsolutions. The optimization was run many times in parallel
until many results produced methods with similar optimality, indicating that we likely obtained values near the true
minimum.
The parameters N and M are the bounds on the stability
region and also represent a trade-off between the stability inthe drift and the stability in the diffusion. A method which
is optimized when M is small would be highly stable inthe case of small noise, but would not be guaranteed tohave good stability properties in the presence of large noise.
Thus these parameters are knobs for tuning the algorithms
for specific situations, and thus we solved the problem fordifferent combinations of N and M to determine different
algorithms for the different cases.

3.1. The Stability Equation for Order 1.5 SRK
Methods with Diagonal Noise
For diagonal noise, we use the mean-square definition of stability[26].
 A method is mean-square stable if
2
limn→∞ E |Xn | = 0 on the test equation
dXt = µXt dt + σXt dWt .

(8)

In matrix form we can re-write our method as given by
Xn+1 = Xn + µh α · H (0) + σI(1) β (1) · H (1)
(9)




I(1,0)
I(1,1)
β (3) · H (1)β (2) · H (1) + σ
+σ √
h
h
(10)

I(1,1,1)  (4)+σβ · H (1)(11)
hwith stages
H (0)
H (1)

I(1,0) (0) (1)B H ,(12)h√= Xn + µ∆tA(1) H (0) + σ ∆tB (1) H (1)= Xn + µ∆tA(0) H (0) + σ

where Xˆn is the size s constant vectorh 2of iXn . Let z = µh√
U
. In this space,
and w = σ h and isolate S = E Un+12nz is the stability variable for the drift term and w is the√ stability in the diffusion term. Under this scaling h, h ,the equation becomes independent of h and thus becomesa function S(z, w) on the coefficients of the SRK method
where mean-square stability is achieved when |S(z, w)| < 1.

3.2. An Optimization Problem for Determinationof Coefficients
We wish to determine the coefficients for the diagonal
SRK methods which optimize the stability. To do so, we
generate an optimization problem which we can numericallysolve for the coefficients. To simplify the problem, we letz, w ∈ R. Define the function
Z M Z 1f (z, w; N, M ) =χS(z,w)≤1 (z, w)dzdw. (13)
−M

−N

max

f (z, w)

A(i) ,B (i) ,β (i) ,α(14)subject to: Order Constraints

3.3. Resulting Approximately-Optimal Methods
The coefficients generated for approximately-optimal
methods fall into three categories. In one category we havethe drift-dominated stability methods where large N andsmall M was optimized. On the other end we have thediffusion-dominated stability methods where large M andsmall N was optimized. Then we have the mixed stability
methods which used some mixed size choices for N and
M . As a baseline, we optimized the objective without
constraints on the ci to see what the “best possible method”
would be. When this was done with large N and M , theresulting method, which we name SOSRI, has almost every(0)value of c satisfy the constraints, but with c2 ≈ −0.04 and(0)c4 ≈ 3.75. To see if we could produce methods which were(B)Euler-Maruyama

w
3

2

2

1
- 10

- 5

- 10

dt

- 5

- 1

-1

- 2

-2

(C)
Figure 3. Convergence results on Equation 15. The error is averaged over1000 trajectories. Shown are the strong l2 error along the time series ofthe solution. The test used a fixed time step h = 1/2−4 to h = 1/2−7 .

-3

(E)

SOSRI

- 10

SOSRI
SOSRI2
1.5 Order Reference

1
z

- 3

(D)

Linear Convergence Tests

SRIW1

w
3

Error

(A)

SOSRI2w
3w
3

4. Numerical Results22

4.1. Convergence Tests

1

1
z

- 5

- 10

- 5

- 1

-1

- 2

-2

- 3

-3

Figure 2. SOSRI Stability Regions. The stability regions(S(z, w) ≤ 1)for the previous and SOSRI methods are plotted inthe (z, w)-plane. (A) Euler-Maruyama. (B) SRIW1. (C) SRIW2. (D)
SOSRI. (E) SOSRI2

more diffusion-stable, we decreased N to optimize more in
w but failed to produce methods with substantially enlarged
diffusion-stability over SOSRI.
Adding only the inequality constraints on the ci andlooking for methods for drift-dominated stability, we failed
to produce methods whose ci estimators adequately covered
the interval. Some of the results did produce stability regions(0)
similar to SOSRI but with ci < 0.5 which indicates themethod could have problems with error estimation. Whenplacing the equality constraints on the edge ci , one method,which we label SOSRI2, resulted in similar stability to
SOSRI but satisfy the ci constraints. In addition, this method
(1)
(1)
(0)
(0)
satisfies c3 = c4 = 1 and c3 = c4 = 1. The stability
regions for these methods is shown in Figure 2.
To look for more diffusion-stable methods, we dropped
to N = 6 to encourage the methods to expand the stability
in the w-plane. However, we could not find a method whose
stability region went substantially beyond [−2, 2] in w. Thiswas further decreased to N = 1 where methods still could
not go substantially beyond |2|. Thus we were not able toobtain methods optimized for the diffusion-dominated case.
This hard barrier was hit under many different constraint
and objective setups and under thousands of optimization
runs, indicating there might be a diffusion-stability barrier
for explicit methods.

In order to test the efficiency and correctness of the SRI
algorithms, we used the linear test Appendix Equation 15.
Figure 3B demonstrates that the SOSRI methods achieve thestrong order 1.5 on Equation 15.

4.2. Epithelial-Mesenchymal Transition
Model (20 Pathwise Stiff SDEs)(EMT)

To test the real consequences of the enhanced stability,
we use the Epithelial-Mesenchymal Transition (EMT) modelof 20 pathwise stiff reaction equations introduced in [3],studied as a numerical test in [9], and written in Appendix
Section B. In the previous work it was noted that t ∈ [0, 1]
was a less stiff version of this model. Thus we first tested the
speed that the methods could solve for 10,000 trajectories
with no failures due to numerical instabilities. The toleranceswere tuned for each method by factors of 2 and finding thelargest values that were stable. Since SOSRI demonstrated
that its stability is much higher than even SOSRI2, we show
the effect of tolerance changes on SOSRI as well. The results
show that at similar tolerances the SOSRI method takes
nearly 5x less time than SRIW1 (Table 1). However, there
is an upper bound on the tolerances before the adaptivity is
no longer able to help keep the method stable. For SRIW1,
this bound is much lower, causing it to run more than 15x
slower than the fastest SOSRI setup. Interestingly SOSRI2
required a higher tolerance than SRIW1 but was 3x fasterthan SRIW1’s fastest setup. We note that SOSRI’s highest
relative tolerance 2−7 ≈ 7 × 10−3 is essentially requiring
4 digits of accuracy (in strong error) when considering
the conservativeness of the error estimator, which is far
beyond the accuracy necessary in many cases. Lastly, we
note that the SOSRI method is able to solve for 10,000 stabletrajectories more than 60x faster than any of the tested fixed
time step methods.
We then timed the run time to solve 10 trajectories inthe t ∈ [0, 500] case (Table 2). This time we found the
optimal tolerance in terms of powers of 10. Once again,

Algorithm

AbstolReltol

Run-time (seconds)

Relative Time (vs SOSRI)

Algorithm

AbstolReltol

Run-time (seconds)

Relative Time (vs SOSRI)

SOSRI

2−7
2−7

2−4
2−6

2.62

1.0x

SOSRI

1.0x

1.0x

2−12
2−13

2−15
2−7

8.78

3.3x

SOSRI2Euler-Maruyama

10−2
10−3

700.76

2.75

10−1
10−3

1016.61
1758.85

1.5x
2.5x

3.05

1.2x

SRIW1

10−5

10−3

4205.52

6.0x

2−12
2−13

2−15
2−11

8.69

3.3x

5.56

2.2x

2−13

2−7

15.16
169.96
182.59424.30

5.8x
64.8x
69.6x
161.7x8912.91

3396x

SOSRI
SOSRI
SOSRI
SOSRI2SOSRI2
SRIW1
Euler-Maruyama
Runge-Kutta Milstein
Fixed Time-step
SRIW1


DISTM θ = 1
2

TABLE 1. SRI TIMES FOR THE THE EMT MODEL ON t ∈ [0, 1]. T HE
EQUATIONS WERE SOLVED 10,000 TIMES WITH THE GIVENTOLERANCES TO COMPLETION AND THE ELAPSED TIME WASRECORDED . T HE FIXED TIME STEP METHODS HAD THEIR ∆t
DETERMINED AS THE LARGEST ∆t IN INCREMENTS OF POWERS OF 2
THAT PRODUCED NO UNSTABLE TRAJECTORIES , AS SHOWN IN [9].
DISTM IS THE D RIFT-I MPLICIT S TOCHASTIC θ -M ETHOD

Algorithm

AbstolReltol

Run-time (seconds)

Relative Time (vs SOSRI)

SOSRI

10−2
10−4

10−2
10−4

22.47

1.0x

73.62

3.3x

10−5
10−4

10−3
10−4

89.19

4.0x

76.12

3.4x

10−5
10−5

10−3
10−3

121.75

5.4x

SOSRI
SOSRI
SOSRI2SOSRI2SRIW1DIRKM θ = 1
2
DIEM θ = 1
2

147.896.6x

7378.55

328.3x8796.47

391.4xTABLE 2. SRI TIMES FOR THE THE EMT MODEL ON t ∈ [0, 500].
T HE EQUATIONS WERE SOLVED 10 TIMES WITH THE GIVENTOLERANCES TO COMPLETION AND THE ELAPSED TIME WASRECORDED . T HE FIXED TIMESTEP METHODS HAD THEIR ∆t CHOSEN
BY INCREMENTING BY 10−5 UNTIL 10 CONSECUTIVE TRAJECTORIESWERE STABLE . D RIFT-I MPLICIT E ULER M ARUYAMA (DIEM) HAD
1
AND D RIFT-I MPLICIT RUNGE -K UTTA M ILSTEIN∆t = 60000
1(DIRKM) HAD ∆t = 50000.

SRIW1 needed a lower tolerance than is necessary in orderto stay stable. SOSRI is able to solve the problem only
asking for around tol = 10−2 , while the others require more(especially in absolute tolerance as there is a stiff reactant
whose values travel close to zero). One interesting point tonote is that at similar tolerances both SOSRI and SOSRI2
receive similar timings and both over 6 times faster than
the fastest SRIW1 tolerance setup. Both are nearly twice
as fast as SRIW1 when matching tolerances as well. Given
the conservativeness of the error estimators generally beingaround 2 orders of magnitude more precise than the local error estimate, the low tolerance solutions are accurate enoughfor many phenomenological experiments and thus present a
good speedup over previous methods. The timings for EulerMaruyama and Runge-Kutta Milstein schemes are omitted
since the tests were unable to finish. From the results of [9]we note that the average dt for SRIW1 on the edge of its
stability had that the smallest dt was approximately 10−11 .
The stability region for fixed step-size Euler-Maruyama isstrictly smaller than SRIW1 (Figure 2) which suggests that it
would require around 5×1012 time steps (with Runge-Kutta
Milstein being similar) to solve to t = 500. Thus, given it
takes on our setup extrapolating the time given 170 seconds
for 220 steps, this projects to around 1.6 × 108 seconds, orapproximately 5 years.

TABLE 3. SRI TIMES FOR THE THE RETINOIC ACID SPDE MODEL
ON t ∈ [0, 500] . T HE EQUATIONS WERE SOLVED TWICE WITH THE
GIVEN TOLERANCES TO COMPLETION AND THE ELAPSED TIME WASRECORDED . T HE TOLERANCES WERE CHOSEN AS THE HIGHEST PAIR OF
TOLERANCES WHICH DID NOT DIVERGE ( GOING UP BY POWERS OF 10).
N OTE THAT NONE OF THE CASES DID THE TWO TIMINGS VARY BY
MORE THAN 1% OF THE TOTAL RUN TIME . E ULER -M ARUYAMA USEDTIME STEPS OF ∆t = 1/20000 SINCE WE NOTE THAT AT ∆t = 1/10000
APPROXIMATELY HALF OF THE TRAJECTORIES ( SIMULATING 10) WEREUNSTABLE .

4.3. Retinoic Acid Stochastic Partial Differential
Equation Model (6x20x100 Semi-Stiff SDEs)
As another test we applied the methods to a methodof lines discretization of a stochastic partial differential
equation (SPDE) describing the spatial regulation of theZrafish hindbrain via retinoic acid signaling ( Section C)[1]. The discretization results in a system of 6 × 20 × 100
SDEs. Starting from an initial zero state, a concentration
gradient emerges over t ∈ [0, 500]. Each of the methods
solved the problem at the highest tolerance that was stablegiving the results in Table 3. Time stepping for this problem
is heavily limited by the high diffusion constant which
results in a strict CFL condition for the 2nd order finite
difference discretization that is used (in the PDE sense),making this problem’s stepping stability-bound for explicit
methods. Because of this stiffness in the real axis, wefound that the previous high order adaptive method SRIW1
did not perform well on this problem in comparison toEuler-Maruyama because the drift term is expensive and theextra function calls outweighed the slightly larger timesteps.
However, the enhanced stability of the SOSRI and SOSRI2
methods allowed for much larger time steps while keeping
the same number of f calls per step, resulting in a more
efficient solution when high accuracy is not necessary. We
note that the drift-implicit stochastic θ-method and drift implicit θ Runge-Kutta Milstein methods were too inefficientto estimate since their time steps were constrained to be nearthat of the Euler-Maruyama equation due to divergence of
the Newton iterations.

5. Discussion
In this work we derived stability-optimized SRK methods for diagonal noise equations. Importantly, our derivation methods utilized heavy computational tools in order
to approximately optimize otherwise intractable equations.
This same method of derivation can easily be scaled up
to higher orders, and by incorporating the coefficients forhigher conditions, efficiency can be optimized as well byadding the norm of the principle error coefficients to the
optimization function. The majority of the search was performed using global optimizers in massive parallel using ahand-optimized CUDA kernel for the numerical integral of
the characteristic function, replacing man-hours with corehours and effectively optimizing the method. The clear next
steps are to find SRI methods with minimal error estimatesand sensible stability regions for the cases in which lower
strong error matters, and similar optimizations on SRK
methods developed for small noise problems. We note thathigh strong order methods were investigated because of
their better trajectory-wise convergence, allowing for a morerobust solution and error estimation since our application totransiently pathwise stiff equations requires such properties.
The main caveat for our methods is the restrictions
on the form of noise. Further research should focus on
the expansion of this these techniques to high order nondiagonal noise integrators. In addition, when g is non-zeroa “diagonal noise” problem over the complex plane does nothave diagonal noise (due to the mixing of real and complex
parts from complex multiplication, and reinterpretation as a2n real system). Thus these methods are not applicable toproblems defined in the complex plane with complex Wiener
processes. Development of similar integrators for commutative noise problems could allow for similar performance
benefits on such problems and is a topic for future research.
Our timings show that the current high order SRK
methods are stability-bound and that when scientific studies
are only looking for small amounts of accuracy in stochastic
simulations, most of the computational effort is lost to generating more accurate than necessary solutions in order to satisfy stability constraints. For diagonal noise approximately
6x faster than the current adaptive methods (SRIW1), whilecommon methods like Euler-Maruyama and Drift-Implicit θ
Runge-Kutta Milstein were in many cases hundreds of times
slower or in many cases could not even finish. We have also
shown that these methods are very robust even at high tolerances and have a tendency to produce the correct qualitative
results on semi-stiff equations (via plots) even when the userchosen accuracy is low. Given that the required user input
is minimal and work over a large range of stiffness, we seethese as very strong candidates for default general purpose
solvers for problem-solving environments such as MATLAB
and Julia since they can easily and efficiently produce results
which are sufficiently correct.

to a chemical reaction network modeled via mass-action
kinetics with Hill functions for the feedbacks. This model
was introduced in [3].
Ad [S1]t−n0
n1
S + [OV OL2] /J1
SS([T GF ] + [T GF 0] (t)) /J0S n0S
k0S + kS
(1 + A) 1 + [S] /J2S
kdS ([S1] − [SR]) − kdSR [SR]=kS ([S1] − [SR]) − kdS [S]=kO34 +=([T GF ] + [T GF 0] (t)) /J0S

=

dt
d [S]
dt
d [miR34]−=

T k KSR ([S1] − [SR]) ([miR34] − [SR]) − [SR]d [SR]dt

Z[S] /J1Z
1 + [S] /J1Z n1Z + [OV OL2] /J2Z n2Z

5
5
X
X

i 
i
kdZR C5 ZRi
C5 [ZR] −kdZ [Z] −
i
i=1
i=1


5
X

i 
C5 ZRi  − kdZ [Z]kZ [Z] −i=1
k200k0200 +
1 + ([S] /J1200 )n1200 + ([Z] /J2200 )n2200


5
X

i 
kd200 [mir] −iC5 ZRi − [T R]
i=1

d [Z]
=

k0Z + kZ

dt



−
d [Z]
=
dt
d [mir]
=
dt
−

5
X




i 
1 − λi kdZR C5 i ZRi − 1 − λT R kdT R [T R]
i
i=1



5
X

i 
T k K1 [mir] −iC5 ZRi − [T R]
i=1



5
X

i 
[Z] −C5 ZRi  − [ZR1 ]
i=1




5
X

i 

T k K2 [mir] −iC5 ZRi − [T R] [ZR1 ] − [ZR2 ]
i=1




5
X

i 
T k K3 [mir] −iC5 ZRi − [T R] [ZR1 ] − [ZR3 ]

−d [ZR1 ]=dt

d [ZR2 ]=dt
d [ZR3 ]

=

dt

i=1

d [ZR4 ]=

ktgf − kdtgf ([tgf ] − [T R]) − kdT R [T R]=

k0T GF + kT GF ([tgf ] − [T R]) − kdT GF [T GF ]




5
X

i 
T k KT R [mir] −iC5 ZRi − [T R] ([tgf ] − [T R]) − [T R]
i=1
kE1
kE2
k0E +
+
− kdE [Ecad]
1 + [S] /J1E n1E
1 + [Z] /J2E n2E
n1
n2
V
V
[S] /J1V[Z] /J2V
kV 1+ kV 2
1 + [S] /J1V n1V
1 + [Z] /J2V n2VBk0V + − kdV [V im]
1 + [OV OL2] /J3V
1
k00 + k0
− kdO [OV OL2]
1 + ([Z] /J0 )n0

i=1



dt
d [T GF ]dt
d [T R]=dt
d [Ecad]
=
dt
=

=

=
dt
d [OV OL2]p

Xt = X0 e.

2. Epithelial-Mesenchymal Transition Model
The Epithelial-Mesenchymal Transition (EMT) model isgiven by the following system of SDEs which correspond

5
X



where

i 
iC5 ZRi=i 
C5 ZRi

=

i=1
5
X




i 
iC5 ZRi − [T R] [ZR1 ] − [ZR5 ]

kOp [OV OL2] − kdOp [OV OL2]p

=

dt

(16)



i=1

d [tgf ]5X

T k K5 [mir] −

dt



2
β− α2 t+αWt




i 
iC5 ZRi − [T R] [ZR1 ] − [ZR4 ]

=d [OV OL2]

(15)

5
X

T k K4 [mir] −dtAppendix=dtd [ZR5 ]B

1dXt = αXt dt + βXt dWt X0 = ,2
11
where α = 10
and β = 20
with true solution

n1d [V im]1. Diagonal Noise Test Equationk34

1 + ([S] /J134 )n134 + ([Z] /J234 )n234
kd34 ([miR34] − [SR]) − 1 − λSR kdSR [SR]dt

5 [ZR1 ] + 20 [ZR2 ] + +30 [ZR3 ] + 20 [ZR4 ] + 5 [ZR5 ] ,5 [ZR1 ] + 10 [ZR2 ] + 10 [ZR3 ] + 5 [ZR4 ] + [ZR5 ] ,i=1([T GF 0] (t)

=

1
2
0

t > 100
o.w.

The parameter values are given in Table 4.

Parameter

Value

Parameter

Value

Parameter

Value

Parameter

Value

J1200

3

J1E

0.1

K2

1

k0O

0.35

J2200

0.2

J2E

0.3

K3

1

kO200

0.0002

J134

0.15

J1V

0.4

K4

1

kO34

0.001

J234

0.35

J2V

0.4

K5

1

kdS

0.09

JO

0.9

J3V

2

KT R

20

kdtgf

0.1

J0S

0.6

J1Z

3.5

KSR

100

kdZ

0.1

J1S

0.5

J2Z

0.9

T GF 00

kdT GF

0.9

J2S

1.8

K1

1

Tk

1000

kdZ

1.66

k0S

0.0005

k0Z

0.003

λ1

0.5

k0T GF

1.1

n1200

3

n1S

2

λ2

0.5

k0E

5

n2200

2

n1E

2

λ3

0.5

k0V

5

n134

2

n2E

2

λ4

0.5

kE1

15

n234

2

n1V

2

λ5

0.5

kE2

5

nO

2

n2V

2

λSR

0.5

kV 1

2

n0S

2

n2Z

6

λT R

0.5

kV 2

5

kO

1.2

k200

0.02

k34

0.01

ktgf

0.05

kZ

0.06

kT GF

1.5

kS

16

kZ

16

kdZR
1
kdZR
5
kdSR

0.5

kdZR
2
kdO

0.5

kdZR

0.5

kdZR
4
kd34

0.5

0.5

1.0

3

kd200

0.035

4. SOSRI
Coefficient
(0)
A2,1
(0)
A3,1
(0)
A3,2
(0)
A4,1
(0)
A4,2
(0)
A4,3
(1)
A2,1
(1)
A3,1
(1)
A3,2
(1)
A4,1
(1)
A4,2
(1)
A4,3
(0)
B2,1
(0)
B3,1
(0)
B3,2
(0)
B4,1
(0)
B4,2
(0)
B4,3
(1)
B2,1
(1)
B3,1
(1)
B3,2
(1)
B4,1
(1)
B4,2
(1)
B4,3

0.035

kOp
10
kdOp
10
TABLE 4. TABLE OF PARAMETER VALUES FOR THE EMT M ODEL .
0.9

Parameter

kdE

0.05

Value

kdV

Parameter

0.05

Value

Parameter

Value

σRAin ,σDN A ,σRAout
0.1
ω
100
u
0.01
b
0.17
γ
3.0
d
0.1
α
10000
δ
0.0013e1
β01
η0.0001a1
c0.1
r0.0001
ζ
0.02
ν
0.85
λ
0.85
D
250.46TABLE 5. TABLE OF PARAMETER VALUES FOR THE RA SPDE M ODEL .

α1α2

3. Retinoic Acid SPDE Model



out
d [RAout ] = β(x) + D∆ [RAout ] − b [RAout ] + c RAin dt + σRA
dWt
out
!
!
α [DN A]d RAin = b [RAout ] + δ [BP ] [DN A] − γ [BP ] + η +
−c
RAin
dt
ω + [DN A]d [RABP ] = γ [BP ] RAin + λ [BP ] [DN A] − (δ + ν [RAR]) [RABP ] dt
DN A
d [DN A] = (ν [RABP ] [RAR] − λ [BP ] [DN A]) dt + σDN A [DN A] dWtd [BP ] = a − λ [BP ] [DN A] − γ [BP ] RAin + (δ + ν [RAR]) [RABP ] dt!
d [DN A]+ −u [BP ] +dt
e + [DN A]

where β(x) = β0 H(x − 40) with H the Heaviside step
function and x = 40 is the edge of retinoic acid production[1]. The space was chosen as [−100, 400] × [0, 100] with∆x = ∆y = 5. The boundary conditions were no-flex on
every side except the right side which had leaky boundary
conditions with parameter kA = 0.002, though full no-flux
does not noticeably change the results. The parameter values
are given in Table 5.
All entries not listed are zero.

Coefficient

Value

α3

0.4736296532772559

2.842612915017106

α4
(0)
c2

0.026404498125060714-2.0527723684000727
4.338237071435815(0)c3

0.7898405466170333-2.8895936137439793(0)c4

3.7504010171562823

2.3017575594644466

(1)
c1

0

0.26204282091330466(1)
c2

0.26204282091330466

0.20903646383505375(1)
c3

0.05879875232001766-0.1502377115150361(1)c4

0.05836595312746999
0.6149440396332373
0.08535117634046772-0.21641093549612528
1.5336352863679572
0.26066223492647056-1.0536037558179159
1.7015284721089472-0.20725685784180017-0.5119011827621657
2.67767339866713-4.9395031322250995
0.15580956238299215
3.2361551006624674-1.4223118283355949
1.140099274172029-0.6401334255743456

(1)
β1
(1)
β2
(1)
β3
(1)
β4
(2)
β1
(2)
β2
(2)
β3
(2)
β4
(3)
β1
(3)
β2
(3)
β3
(3)
β4
(4)
β1
(4)
β2
(4)
β3
(4)
β4

-0.04199224421316468

0.758661169101175-1.8453464565104432
2.688764531100726-0.2523866501071323
0.40896857551684956
0.4969658141589478-0.5771202869753592-0.12919702470322217
0.2093514975196336
2.8453464565104425-2.688764531100725
0.2523866501071322-0.40896857551684945
0.11522663875443433-0.57877086147738
0.2857851028163886
0.17775911990655704

5. SOSRI2
Coefficientd [RAR] = (ζ − ν [RABP ] [RAR] + λ [BP ] [DN A] − r [RAR]) dt

Value-0.04199224421316468

(0)
A2,1
(0)
A3,1
(0)
A3,2
(0)
A4,1
(0)
A4,2
(0)
A4,3
(1)
A2,1
(1)
A3,1
(1)
A3,2
(1)
A4,1
(1)
A4,2
(1)
A4,3
(0)
B2,1
(0)
B3,1
(0)
B3,2
(0)
B4,1
(0)
B4,2
(0)
B4,3
(1)
B2,1
(1)
B3,1
(1)
B3,2
(1)
B4,1
(1)
B4,2
(1)
B4,3
α1
α2

Value

CoefficientValue

0.13804532298278663

α30.686995463807979

0.5818361298250374
0.4181638701749618

α4
(0)
c2

0.13804532298278663

0.4670018408674211

(0)
c3

1

0.8046204792187386

(0)
c4

1

-0.27162232008616016

(1)
c1

0

0.45605532163856893(1)
c2

0.45605532163856893

0.7555807846451692

(1)
c3

1

0.24441921535482677(1)c4

0.6981181143266059
0.3453277086024727-0.04344582292908241
0.08852381537667678
1.0317752458971061
0.4563552922077882
1.73078280444124-0.46089678470929774-0.9637509618944188
0.6753186815412179-0.07452812525785148-0.49783736486149366-0.5591906709928903
0.022696571806569924-0.8984927888368557-0.15036858140642623
0.7545275856696072

(1)
β1
(1)
β2
(1)
β3
(1)
β4
(2)
β1
(2)
β2
(2)
β3
(2)
β4
(3)
β1
(3)
β2
(3)
β3
(3)
β4
(4)
β1
(4)
β2
(4)
β3
(4)
β4

-0.2911544680711602

1-0.45315689727309133
0.8330937231303951
0.3792843195533544
0.24077885458934192-0.4994383733810986
0.9181786186154077-0.25613778661003145-0.16260245862427797
1.4531568972730915-0.8330937231303933-0.3792843195533583-0.24077885458934023-0.4976090683622265
0.9148155835648892-1.4102107084476505
0.9930041932449877

Acknowledgments
We would like to thank the members of SciML, specifically David Widmann (@devmotion), Yingbo Ma (@Ying-boMa) and (@dextorious) for their contributions to theecosystem. Their efforts have helped make the developmentof efficient implementations possible. This work was partly
supported by the NSF grant DMS1763272 and a grant fromthe Simons Foundation (594598, QN). This work used theExtreme Science and Engineering Discovery Environment(XSEDE), which is supported by National Science Foundation grant number ACI-1548562.

References[16] T. Schaffter, “From genes to organisms: Bioinformatics system models and software,” Thesis, École polytechnique fédérale de Lausanne,
2014.
[17] H. Gilsing and T. Shardlow, “Sdelab: A package forsolvingstochastic
differential
equationsinmatlab,”Journal of Computational and Applied Mathematics, vol.
205, no. 2, pp. 1002–1018, 2007. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0377042706004195[18] A. Janicki, A. Izydorczyk, and P. Gradalski, Computer Simulationof Stochastic Models with SDE-Solver Software Package. Berlin,
Heidelberg: Springer Berlin Heidelberg, 2003, pp. 361–370.

[1]

C. Rackauckas and Q. Nie, “Mean-independent noise control of cellfates via intermediate states,” iScience, vol. Accepted, 2018.

[19] A. Rossler, “Runge kutta methods for the strong approximationof solutions of stochastic differential equations,” SIAM Journal on
Numerical Analysis, vol. 48, no. 3, pp. 922–952, 2010.

[2]

B. Wang and Q. Qi, “Modeling the lake eutrophication
stochastic ecosystem and the research of its stability,”
Mathematical
Biosciences,
2018.
[Online].
Available:
https://www.sciencedirect.com/science/article/pii/S0025556418301780[20] M. Wiktorsson, “Joint characteristic function and simultaneous
simulation of iterated ito integrals for multiple independent brownian
motions,” The Annals of Applied Probability, pp. 470–487, 2001.
[Online]. Available: http://projecteuclid.org/euclid.aoap/1015345301[3]

T. Hong, K. Watanabe, C. H. Ta, A. Villarreal-Ponce, Q. Nie, and
X. Dai, “An ovol2-zeb1 mutual inhibitory circuit governs bidirectionaland multi-step transition between epithelial and mesenchymal states,”
PLoS Comput Biol, vol. 11, no. 11, p. e1004569, 2015.

[21] F. S. Lawrence, “Some practical runge-kutta formulas,” Math. Comput., vol. 46, no. 173, pp. 135–150, 1986.

[4]

[5]

[6]

D. K. Wells, W. L. Kath, and A. E. Motter, “Control ofstochastic and induced switching in biophysical networks,” Phys.
Rev. X, vol. 5, p. 031036, Sep 2015. [Online]. Available:
https://link.aps.org/doi/10.1103/PhysRevX.5.031036
A. K. Ghosh, F. Hussain, S. Jha, C. J. Langmead, and S. K. Jha,
“Discovering rare behaviours in stochastic differential equations using
decision procedures: applications to a minimal cell cycle model,”International Journal of Bioinformatics Research and Applications
2, vol. 10, no. 4-5, pp. 540–558, 2014.
S. Donnet and A. Samson, “A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models,”
Advanced drug delivery reviews, vol. 65, no. 7, pp. 929–939, 2013.

[7]

E. Hairer and G. Wanner, Solving Ordinary Differential Equations II- Stiff and Differential-Algebraic Problems. Springer, 1991.

[8]

T. Li, A. Abdulle, and W. E, “Effectiveness of implicit methods forstiff stochastic differential equations,” Commun. Comput. Phys, vol. 3,
no. 2, pp. 295–307, 2008.

[9]

C. Rackauckas and Q. Nie, “Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling withmemory,” Discrete and Continuous Dynamical Systems - Series B,
vol. 22, no. 7, pp. 2731–2761, 2016.

[10] J. Lawson, “An order five runge-kutta process with extended regionof stability,” SIAM Journal on Numerical Analysis, vol. 3, no. 4, pp.
593–597, 1966.
[11] D. I. K. Ahmadia and A. J., “Optimal stability polynomials fornumerical integration of initial value problems,” CAMCOS, vol. 7,
no. 2, pp. 247–271, 2012.
[12] P. J. van der Houwen, “Explicit runge-kutta formulaswith increased stability boundaries,” Numerische Mathematik,
vol. 20, no. 2, pp. 149–164, 1972. [Online]. Available:
http://dx.doi.org/10.1007/BF01404404[13] Y. Komori and K. Burrage, “Strong first order s-rock methods forstochastic differential equations,” Journal of Computational and Applied Mathematics, vol. 242, no. Supplement C, pp. 261–274, 2013.
[14] A. Abdulle and S. Cirilli, “S-rock: Chebyshev methods forstiff stochastic differential equations,” SIAM Journal on Scientific
Computing, vol. 30, no. 2, pp. 997–1014, 2008. [Online]. Available:
https://doi.org/10.1137/070679375[15] Y. Komori and K. Burrage, “Weak second order s-rock methods forstratonovich stochastic differential equations,” Journal of Computational and Applied Mathematics, vol. 236, no. 11, pp. 2895–2908,
2012.

[22] J. R. Dormand and P. J. Prince, “A family of embedded
runge-kutta formulae,” Journal of Computational and Applied
Mathematics, vol. 6, no. 1, pp. 19–26, 1980. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/0771050X80900133[23] J. C. Butcher, “Numerical methods for ordinary differential equationsin the 20th century,” Journal of Computational and Applied
Mathematics, vol. 125, no. 1-2, pp. 1–29, 2000. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0377042700004556[24] C. Tsitouras, “Runge-kutta pairs of order 5(4) satisfying only the firstcolumn simplifying assumption,” Computers & Mathematics with
Applications, vol. 62, no. 2, pp. 770–775, 2011. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0898122111004706[25] W. H. Enright, D. J. Higham, B. Owren, and P. W. Sharp, “A surveyof the explicit runge-kutta method,” 1995.
[26] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations. Springer Berlin Heidelberg, 2011.
[27] I. Dunning, J. Huchette, and M. Lubin, “Jump: A modeling language
for mathematical optimization,” SIAM Review, vol. 59, no. 2, pp. 295–
320, 2017. [Online]. Available: https://doi.org/10.1137/15M1020575[28] S. G. Johnson, “The nlopt nonlinear-optimization package,” 2018.
[Online]. Available: http://ab-initio.mit.edu/nlopt[29] A. Conn, N. Gould, and P. Toint, “A globally convergent
augmented lagrangian algorithm for optimization with general
constraints and simple bounds,” SIAM Journal on Numerical
Analysis, vol. 28, no. 2, pp. 545–572, 1991. [Online]. Available:
https://doi.org/10.1137/0728030[30] E. G. Birgin and J. M. Martinez, “Improving ultimate convergenceof an augmented lagrangian method,” Optimization Methods andSoftware, vol. 23, no. 2, pp. 177–195, 2008. [Online]. Available:
https://doi.org/10.1080/10556780701577730[31] T. P. Runarsson and Y. Xin, “Search biases in constrained evolutionary
optimization,” IEEE Transactions on Systems, Man, and Cybernetics,
Part C (Applications and Reviews), vol. 35, no. 2, pp. 233–243, 2005.