Appendix C: Numerical solutions of ODES

---

import numpy as np
import scipy.integrate
import scipy.optimize

import bokeh.io
import bokeh.plotting

bokeh.io.output_notebook()

Loading BokehJS ...

---

In Technical Appendix 2a, we learned basic syntax of how to use scipy.integrate.odeint() to perform numerical integration of ordinary differential equations. We gave a cursory description of how it works, but the integration was more or less treated like a black box. Here, we give a more pedagogical introduction to how numerical solution of systems of ODEs works and why it is necessary to use a package like scipy.integrate.odeint().

In our examples, we mostly use dynamical equations for chemical kinetics. Many of the dynamical equations we encounter in the study of biological circuits are equations describing chemical kinetics or are otherwise similar with similar issues.

Euler’s method as an introduction to numerical methods

A classic introduction to numerical methods for solving system of time-based ODEs involves Euler’s method. Consider an ordinary differential equation

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} = f(y, t), \end{align}

with initial condition \(y(0) = y_0\). Our goal is to find \(y(t)\). For some functions \(f(y, t)\), and analytical solution can be found. For example, if \(f = -k y\), we have exponential decay and \(y(t) = y_0\mathrm{e}^{-kt}\). But we cannot always find an analytical solution.

Our goal instead is to try to find \(y(t)\) at discrete time points. For Euler’s method in particular, we will chose time points \(t = \{0, \Delta t, 2\Delta t, \ldots n\Delta t\}\), where \(\Delta t\) is some small time step.

Recalling your first forays into calculus, recall one of the definitions of a derivative,

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} = \lim_{\Delta t\to0}\frac{y(t+\Delta t) - y(t)}{\Delta t}. \end{align}

We will drop the limit part, but keep \(\Delta t\) small, giving us

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} \approx \frac{y(t+\Delta t) - y(t)}{\Delta t}. \end{align}

This expression for the time derivative of \(y\) is referred to as the forward differencing formula. Inserting this into \(\mathrm{d}y/\mathrm{d}t = f(y, t)\) and rearranging, we have

\begin{align} y(t+\Delta t) \approx y(t) + \Delta t f(y(t), t). \end{align}

This is a nice formula, because if we start at time \(t = 0\), we know \(y(0) = y_0\), so we can solve for \(y(\Delta t)\) as

\begin{align} y(\Delta t) \approx y_0 + \Delta t f(y_0, 0). \end{align}

But now that we have \(y(\Delta t)\), we can use the formula again to get \(y(2\Delta t)\).

\begin{align} y(2\Delta t) \approx y(\Delta t) + \Delta t f(y(\Delta t), \Delta t). \end{align}

And we can keep marching forward in time until we get \(y\) at all of the desired time points. This technique is more or less how initial value problems are solved. Starting from an initial condition, we use a differencing formula to figure out what the value of the solution to the ODE(s) is at a time point a little bit in the future. We then keep marching forward in time.

Euler’s method is an example of an explicit technique for solving ODEs, since the value of the solution at time \(t + \Delta t\) depends only on the right-hand side evalued at time \(t\). That is, \(y(t + \Delta t)\) depends only on \(f(y, t)\). We will contrast this with implicit methods later.

Now, we can code up Euler’s method. We will assume vectorial \(y\), meaning that we are solving a system of ODEs. There are no additional considerations for a vectorial \(y\). In the function, note that the keyword argument args is a tuple that gives extra parameters that are passed into the function f.

def solve_euler(f, y0, t_max, delta_t, args=()):
    """Solve an ODE defined by y'(t) = f(y, t, *args)
    using Euler's method."""
    # Initialize output array
    y = np.empty((int(t_max / delta_t) + 1, len(y0)))

    # Insert initial condition
    y[0] = y0

    # Start at t = 0, with index 0
    t = 0.0
    i = 0

    # March forward in time
    while i < len(y) - 1:
        y[i+1] = y[i] + delta_t * f(y[i], t, *args)
        i += 1
        t += delta_t

    return y

An example problem: reversible dimerization

To have a concrete example problem, let us consider the dimerization of hexokinase. Hexokinases phosphorylate six-carbon sugars and are important metabolic enzymes. They exist in yeast cells in high concentration, being present in both monomeric and dimeric forms. The two forms bind glucose differently, so the relative concentrations of the monomeric and dimeric forms are of interest. Denoting A as a hexokinase monomer, the chemical reaction is \(\mathrm{A} + \mathrm{A} \rightleftharpoons AA\), with the forward reaction having rate constant \(k_+\) and the reverse reaction having rate constant \(k_-\). The dynamical equations for this system are

\begin{align} &\frac{\mathrm{d}c_\mathrm{AA}}{\mathrm{d}t} = k_+ c_\mathrm{A}^2 - k_-\,c_\mathrm{AA},\\[1em] &\frac{\mathrm{d}c_\mathrm{A}}{\mathrm{d}t} = -2k_+ c_\mathrm{A}^2 + 2k_-\,c_\mathrm{AA}. \end{align}

We can employ the conservation of mass relation

\begin{align} c_\mathrm{A}^0 = c_\mathrm{A} + 2c_\mathrm{AA}, \end{align}

where \(c_\mathrm{A}^0\) is the total concentration of hexokinase molecules, to reduce the dynamical system to a single equation,

\begin{align} \frac{\mathrm{d}c_\mathrm{A}}{\mathrm{d}t} = -2k_+\,c_\mathrm{A}^2 - k_-\,c_\mathrm{A} + k_-\,c_\mathrm{A}^0. \end{align}

To enable solving this using Euler’s method, we need to code up the expression on the right hand side.

def dimerization_rhs(cA, t, k_plus, k_minus, cA_total):
    """
    Right hand side of dimerization ODE
    """
    return -(k_minus + 2 * k_plus * cA) * cA + k_minus * cA_total

Note that I have explicitly put in time dependence in the function definition although \(t\) does not appear in the expression. I do this because the Euler stepper we wrote above expects this, as does the automated integrator we will use later in this appendix.

To solve, we need to specify parameters. We will use parameters obtained by Hoggett and Kellett, 1992, which they obtained for the PI hexokinase in yeast in a solution containing 10 mM glucose. They determined that \(k_+ = 0.25\) (µM·s)⁻¹ and \(k_- = 3.1\) s⁻¹. In one of their experiments, they used \(c_\mathrm{A}^0 = 10\) µM.

# Parameters
k_plus = 0.25  # 1 / µM-s
k_minus = 3.1  # 1 / s

# Total concentration of hexokinase molecules
cA_total = 10  # µM

# Initial concentration of monomer (assume all monomeric)
cA0 = 10  # µM

# Total time we want to use
t_max = 1  # s

# Time step for integrator (50 ms)
delta_t = 0.05  # s

# Assemble parameters into args
args = (k_plus, k_minus, cA_total)

Now we can solve using Euler’s method!

# Generate time points for plotting
t_coarse = np.arange(int(t_max / delta_t) + 1) * delta_t

# Solve for monomeric concentration
cA = solve_euler(
    dimerization_rhs, np.array([cA0]), t_max, delta_t, args=(k_plus, k_minus, cA_total)
)

# Compute dimeric concentration using conservation of mass
cAA = (cA_total - cA[:, 0]) / 2

Let’s take a look at the plot of our solution.

p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=300,
    x_axis_label="time (s)",
    y_axis_label="concentration (µM)",
    x_range=[0, 1],
)
p.line(t_coarse, cA[:, 0], line_width=2, legend_label="monomer")
p.line(t_coarse, cAA, color="orange", line_width=2, legend_label="dimer")

bokeh.io.show(p)

This seems to capture what we would expect, a decay of the monomer concentration with a corresponding increase in the dimer concentration.

As is often the case when we are testing out numerical methods, it is good to compare to a known, analytical result. In this case, we can calculate the rather complicated analytical result. (We will not show the calculation here; only state the result.) Defining

\begin{align} \zeta = \frac{1}{\sqrt{1 + 8 \,\frac{k_+}{k_-}\,c_\mathrm{A}^0}}, \end{align}

the solution is

\begin{align} c_\mathrm{A}(t) = \frac{1}{4}\frac{k_-}{k_+}\left(\zeta^{-1} \tanh\left[\frac{k_- t}{2\zeta} + \mathrm{arctanh}\left(\zeta\left(1 + 4\, \frac{k_+}{k_-}\,c_\mathrm{A,0}\right)\right)\right] - 1\right) \end{align}

if \begin{align} \zeta\left(1+4\,\frac{k_+}{k_-}\,c_\mathrm{A,0}\right) < 1 \end{align}

and

\begin{align} c_\mathrm{A}(t) = \frac{1}{4}\frac{k_-}{k_+}\left(\zeta^{-1} \coth\left[\frac{k_- t}{2\zeta} + \mathrm{arccoth}\left(\zeta\left(1 + 4\, \frac{k_+}{k_-}\,c_\mathrm{A,0}\right)\right)\right] - 1\right) \end{align}

if

\begin{align} \zeta\left(1+4\,\frac{k_+}{k_-}\,c_\mathrm{A,0}\right) \ge 1. \end{align}

We can code this up in Python, noting that Numpy does not have built-in coth and arccoth functions, so we have write those ourselves.

def coth(x):
    """Hyperbolic cotangent function."""
    return 1 / np.tanh(x)


def arccoth(x):
    """Hyberbolic arccotangent function."""
    return np.log((x + 1) / (x - 1)) / 2


def cA_theor(t, k_plus, k_minus, cA_total, cA0):
    """Analytical result for concentration of monomer over time
    for reversible dimerization."""
    zeta = 1 / np.sqrt(1 + 8 * k_plus / k_minus * cA_total)
    arc_arg = zeta * (1 + 4 * k_plus / k_minus * cA0)

    # Which part of the solution we live on
    fun = np.tanh if arc_arg < 1 else coth
    arcfun = np.arctanh if arc_arg < 1 else arccoth

    prefactor = k_minus / k_plus / 4
    return prefactor * (fun(k_minus * t / 2 / zeta + arcfun(arc_arg)) / zeta - 1)

Now we can compute correct values for \(c_\mathrm{A}\) and overlay them on the plot of the numerical result. We will use open circles for the true values.

cA_analytical = cA_theor(t_coarse, k_plus, k_minus, cA_total, cA0)

p.circle(t_coarse, cA_analytical, fill_color=None)
bokeh.io.show(p)

We see that the numerical solution is missing the true values. We could try to remedy this situation by using a smaller \(\Delta t\).

# Smaller Δt
delta_t = 0.01

# Solve for monomeric concentration
cA = solve_euler(
    dimerization_rhs, np.array([cA0]), t_max, delta_t, args=(k_plus, k_minus, cA_total)
)

# Compute dimeric concentration using conservation of mass
cAA = (cA_total - cA[:, 0]) / 2

# Generate new time points for plotting numerical solution
t = np.arange(int(t_max / delta_t) + 1) * delta_t

# Set up plot and populate with analytical solution
p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=300,
    x_axis_label="time (s)",
    y_axis_label="concentration (µM)",
    x_range=[0, 1],
)

# Populate plot and show
p.line(t, cA[:, 0], line_width=2, legend_label="monomer")
p.line(t, cAA, color="orange", line_width=2, legend_label="dimer")
p.circle(t_coarse, cA_analytical, fill_color=None)

bokeh.io.show(p)

Much better!

Variable step sizes can help speed the solution

If we look at the above solution, we are still not quite hitting the analytical solution at short times, though we do well as the system approaches steady state. It makes sense, then, to take shorter time steps when the solution is changing rapidly and longer time steps when it is changing slowly. We would therefore employ a variable step size algorithm.

Such an algorithm is employed in scipy.optimize.odeint(). It uses the Hindmarsh-Petzold algorithm, as implemented in LSODA as part of ODEPACK. (That word salad is mostly so you know what to enter into search engines to find more information about the algorithm.) We will discuss some other features of the algorithm momentarily.

The interface for scipy.integrate.odeint() is the same as for the Euler integration function we just developed, except in place of delta_t and t_max, you specify for which time points you want the solution evaluated. Because scipy.integrate.odeint() does automatic variable time stepping, it chooses its own \(\Delta t\)’s for each time step. It then interpolates its solutions to give you the values at the time points you ask for.

Let’s put this to use. We will use the same time points we used with Euler integration.

# Integrate ODES
cA = scipy.integrate.odeint(
    dimerization_rhs, np.array([cA0]), t, args=(k_plus, k_minus, cA_total)
)

# Add this one to the plot in red
p.line(t, cA[:, 0], line_width=2, color="tomato")

bokeh.io.show(p)

The scipy.integrate.odeint() function much more closely matches the analytical solution at all time points.

Stiffness presents additional challenges

Systems of ODEs that have a wide spread in time scales of dynamics are said to be stiff. Stiff ODEs offer challenges beyond the accuracy and time-stepping issues we have already explored with Euler’s method. Stiff ODEs have different dynamical modes that occur at vastly different time scales. This usually arises from big differences in coefficients in the dynamical equations.

We will consider a classic example of a stiff chemical reaction system proposed by H. H. Robertson, known as the Robertson problem. The chemical reactions and their rate constants are

\begin{align} &\mathrm{Y} + \mathrm{Z} \longrightarrow \mathrm{X};\;\;\;k_1 = 10^4\text{ M}^{-1}\cdot\text{s}^{-1}, \\[1em] &\mathrm{X} \longrightarrow \mathrm{Y};\;\;\;k_2 = 0.04\text{ s}^{-1}, \\[1em] &\mathrm{Y} + \mathrm{Y} \longrightarrow \mathrm{Z};\;\;\;k_3 = 3\times10^7\text{ M}^{-1}\cdot\text{s}^{-1}. \end{align}

Letting lowercase letters be the concentrations of the corresponding uppercase symbols, the dynamical equations are

\begin{align} &\frac{\mathrm{d}x}{\mathrm{d}t} = k_1\,y\,z - k_2\,x,\\[1em] &\frac{\mathrm{d}y}{\mathrm{d}t} = -k_1\,y\,z + k_2\,x - k_3\,y^2,\\[1em] &\frac{\mathrm{d}z}{\mathrm{d}t} = k_3\,y^2. \end{align}

We see a huge range of rate constant values, ranging from \(10^{-2}\) to \(10^{7}\), a sign of stiffness. This is a nice system to study because we know two facts about it. First, an interesting feature of this system is that

\begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} + \frac{\mathrm{d}y}{\mathrm{d}t} + \frac{\mathrm{d}z}{\mathrm{d}t} = 0, \end{align}

so that \(x + y + z\) is always constant. Second, the steady state, the values of \(x\), \(y\), and \(z\) for which all time derivatives vanish, is unique with \(x = y = 0\), as given by setting \(\mathrm{d}z/\mathrm{d}t = \mathrm{d}x/\mathrm{d}t = 0\). Because \(x + y + z\) is constant, the steady state has \(z\) equal to that constant. Knowing these facts will help us analyze correctness of numerical solutions.

Let us now attempt to solve this problem with Euler’s method, starting from an initial condition of \(x = 1\) M and \(y = z = 0\). With our initial condition, \(x + y + z = 1\)M for all time.

def robertson_rhs(xyz, t, k1, k2, k3):
    """Right-hand-side of Robertson's stiff ODEs."""
    # Unpack variables
    x, y, z = xyz

    dx_dt = k1 * y * z - k2 * x
    dz_dt = k3 * y ** 2
    dy_dt = -dx_dt - dz_dt

    return np.array([dx_dt, dy_dt, dz_dt])


# Time steps
delta_t = 0.001
t_max = 10

# Paramters
k1 = 1e4
k2 = 0.04
k3 = 3e7

# Initial condition
xyz0 = np.array([1, 0, 0])
xyz = solve_euler(robertson_rhs, xyz0, t_max, delta_t, args=(k1, k2, k3))

/var/folders/j_/c5r9ch0913v3h1w4bdwzm0lh0000gn/T/ipykernel_80041/1963278859.py:6: RuntimeWarning: overflow encountered in double_scalars
  dx_dt = k1 * y * z - k2 * x
/var/folders/j_/c5r9ch0913v3h1w4bdwzm0lh0000gn/T/ipykernel_80041/1963278859.py:7: RuntimeWarning: overflow encountered in double_scalars
  dz_dt = k3 * y ** 2
/var/folders/j_/c5r9ch0913v3h1w4bdwzm0lh0000gn/T/ipykernel_80041/1963278859.py:8: RuntimeWarning: invalid value encountered in double_scalars
  dy_dt = -dx_dt - dz_dt

We got a warning that we had overflow! This must mean the solver has made grievous errors, since we must always have \(x + y + z = 1\). We can plot the solution to see what was going on. We will put \(x\) in blue, \(y\) in orange, and \(z\) in red.

# Generate time points for plotting numerical solution
t = np.arange(int(t_max / delta_t) + 1) * delta_t

# Make plot
p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=300,
    x_axis_label="time (s)",
    y_axis_label="conc (M)",
    y_axis_type="log",
    x_range=[0, 10],
    y_range=[1e-9, 1e2],
)

# Thin the points so the file of the plot is not so large
p.line(t[::10], xyz[::10, 0], line_width=2)
p.line(t[::10], xyz[::10, 1], line_width=2, color="orange")
p.line(t[::10], xyz[::10, 2], line_width=2, color="tomato")

bokeh.io.show(p)

The concentration axis is on a log scale, so where the \(y\) (orange) curve is missing, we are getting negative values. Importantly, looking right around 9.34 seconds, we see the solution explode. This is an example of instability that can plague explicit methods with stiff systems.

Implicit methods can handle stiffness, but are slower and harder to implement

In contrast to explicit methods, implicit methods are much less susceptible to the types of catastrophic failures due to stiffness we see with explicit methods. Implicit methods employ a different kind of differencing formula. Recall that for Euler’s method, we approximated the time derivative using the forward differencing as

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} \approx \frac{y(t+\Delta t) - y(t)}{\Delta t}, \end{align}

where \(y\) here is a generic solution to a system of ODEs, not to be confused with the concentration of species Y in the Robertson problem. This approximation of the derivative led to updating the solution at time \(t + \Delta t\) as

\begin{align} \frac{y(t+\Delta t) - y(t)}{\Delta t} \approx f(y(t), t). \end{align}

In the backwards Euler method, we instead use backward differencing to obtain a solution at time \(t + \Delta t\) as

\begin{align} \frac{y(t+\Delta t) - y(t)}{\Delta t} \approx f(y(t+\Delta t), t+\Delta t). \end{align}

Our solution is then updated according to

\begin{align} y(t+\Delta t) = y(t) + \Delta t f(y(t+\Delta t), t+\Delta t). \end{align}

This is now an implicit expression for \(y(t+\Delta t)\), since \(y(t+\Delta t)\) appears (in a possibly nonlinear way) on the right hand side. Because implicit methods require solving a set of nonlinear algebraic equations to update the solution, each step takes considerably more time to evaluate than for an explicit method. Also because of this requirement to solve algebraic systems at each time step, they are more difficult to implement. The (very big) upside is that they handle stiffness much better. (We will not go into details as to why and how they perform with respective to stiffness and stability, but that is a big, beautiful field of study in numerics.)

The solution of the above algebraic equation for \(y(t+\Delta t)\) is a from of a fixed point problem of the form

\begin{align} x = g(x). \end{align}

The scipy.optimize module offered a function scipy.optimize.fixed_point() (docs) to solve this problem. We can use this function to code up a backward Euler method.

def solve_backward_euler(f, y0, t_max, delta_t, args=()):
    """Solve an ODE defined by y'(t) = f(y, t, *args)
    using backward Euler's method."""
    # Initialize output array
    y = np.empty((int(t_max / delta_t) + 1, len(y0)))

    # Insert initial condition
    y[0] = y0

    # Start at t = 0, with index 0
    t = 0.0
    i = 0

    # Function to be used in fixed point iteration
    fp_func = lambda y, y_prev, t, delta_t, args: y_prev + delta_t * f(y, t, *args)

    # March forward in time
    while i < len(y) - 1:
        # Package args for fixed point function
        fp_args = (y[i], t, delta_t, args)

        # Solve for new y value and iterate
        y[i + 1] = scipy.optimize.fixed_point(fp_func, y[i], args=fp_args)
        i += 1
        t += delta_t

    return y

Now, let’s try the same calculation as before with backward Euler.

# Solve using backward Euler
xyz = solve_backward_euler(robertson_rhs, xyz0, t_max, delta_t, args=(k1, k2, k3))

# Make same plot as we did with explicit method
p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=300,
    x_axis_label="time (s)",
    y_axis_label="conc (M)",
    y_axis_type="log",
    x_range=[0, 10],
    y_range=[1e-9, 1e2],
)
p.line(t[::10], xyz[::10, 0], line_width=2)
p.line(t[::10], xyz[::10, 1], line_width=2, color="orange")
p.line(t[::10], xyz[::10, 2], line_width=2, color="tomato")

bokeh.io.show(p)

The numerical stability issues are gone. The cost, though, is that the calculation takes much longer per step. We only integrated for 10 seconds, but the level of \(z\) is still very far from its steady state value, so we have a lot more integrating to do! Our implementation of backward Euler would take too long to integrate.

The Hindmarsh-Petzold algorithm is the best of both worlds

As we have seen some of the challenges of numerically solving systems of ODEs, we desire and ODE solver that can:

1. Take small steps when the dynamics are changing rapidly and large steps when they are not.

2. Use explicit time stepping when stiffness is not an issue (for speed), and implicit when it is (for stability).

The Hindmarsh-Petzold algorithm does this. It uses variable time stepping and also detects stiffness and switches to implicit time stepping when necessary.

Importantly, in the Robertson system, we want to take tiiiiiny time steps at the beginning to accurately capture fast transients. These are missing in our naive Euler and backwards Euler time stepping, but scipy.integrate.odeint() captures them. Also because of the variable time stepping and speed, we can run the solver for much much longer, capturing the dynamics all the way up to the steady state. This is really powerful, as it is quite often the case that the dynamics of chemical reaction systems vary across many time scales.

# Time points we want, on a log scale
t = np.logspace(-5, 7, 4000)

# Integrate ODES
xyz = scipy.integrate.odeint(robertson_rhs, xyz0, t, args=(k1, k2, k3))

p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=300,
    x_axis_label="time",
    y_axis_label="conc",
    x_axis_type="log",
    y_axis_type="log",
    y_range=[1e-9, 1e2],
    x_range=[1e-5, 1e7],
)
p.line(t, xyz[:, 0], line_width=2)
p.line(t, xyz[:, 1], line_width=2, color="orange")
p.line(t, xyz[:, 2], line_width=2, color="tomato")

bokeh.io.show(p)

This is great! We capture the short timescale rise of \(y\) (the orange curve) on the scale of tens or milliseconds, and we still see the approach to steady state, which happens on the scale of \(10^5\) seconds, or about one day.

The small kink in the \(z\)-curve (the red one) is from when the solver shifted between explicit to implicit. The error therein is small, though, less than one part per million of the steady state level.

So, going forward, scipy.integrate.odeint() is a very powerful tool for numerically solving chemical systems. For situations where it has troubles, you would need to dig much deeper into the structure of the problem to choose or develop the appropriate numerical methods.

Afterthought: A note about the solving the backward difference equation

Generally, the methods for solving \(y(t+\Delta t) = y(t) + \Delta t f(y(t+\Delta t), t+\Delta t)\) involve techniques in numerical optimization, classically including fixed point iterations, but for illustrative purposes, I show how to solve for time steps with backward Euler for the Robertson system.

For ease of notation, temporarily let \(x = x(t+\Delta t)\), \(y = y(t+\Delta t)\), \(z = z(t+\Delta t)\), and \(x_0 = x(t)\), \(y_0 = y(t)\), \(z_0 = z(t)\). Then, to update the solution from \((x_0, y_0, z_0)\) tp \((x, y, z)\), we have to solve

\begin{align} &x = x_0 + \Delta t (k_1 yz - k_2 x),\\[1em] &y = y_0 + \Delta t(-k_1 yz + k_2 x - k_3 y^2),\\[1em] &z = z_0 + \Delta t k_3 y^2. \end{align}

Inserting the expression for \(z\) into that for \(x\) yields

\begin{align} x = x_0 + \Delta t(k_1 y(z_0 + \Delta t k_3 y^2) - k_2 x). \end{align}

This can be rearranged to give \(x\) as a function of \(y\).

\begin{align} x = \frac{1}{1 + \Delta t k_2}\left(x_0 + \Delta t k_1 y(z_0 + \Delta t k_3 y^2)\right) \end{align}

Inserting this expressions we have derived for \(x\) and \(z\) as functions only of \(y\) into our expression for \(y\) yields

\begin{align} y = y_0 + \Delta t\left(-k_1 y(z_0 + \Delta t k_3 y^2) + \frac{k_2}{1 + \Delta t k_2}\left(x_0 + \Delta t k_1 y(z_0 + \Delta t k_3 y^2)\right) - k_3 y^2\right). \end{align}

This is a cubic polynomial in \(y\) that can be rearranged to give

\begin{align} \Delta t^2\,k_1\,k_3\left(1 - \frac{\Delta t k_2}{1 + \Delta t k_2}\right)y^3 + \Delta t k_3 y^2 + \left(1 + \Delta t k_1 z_0\left(1 - \frac{\Delta t k_2}{1 + \Delta t k_2}\right)\right)y - y_0 -\frac{\Delta t k_2}{1 + \Delta t k_2}\,x_0 = 0. \end{align}

By the Decartes Rule of Signs, this cubic has a single positive real root, which we can find for \(y\), and then determine \(z\) and \(y\) by substitution of the solution for \(y\). This is implemented below.

def solve_timestep_robertson(xyz0, delta_t, k1, k2, k3):
    """Solves for the implicit time step for use in backward
    Euler for the Robertson problem."""
    # Unpack previous time step
    x0, y0, z0 = xyz0

    # Compute coefficients of cubic
    k2_ratio_term = delta_t * k2 / (1 + delta_t * k2)
    a = delta_t**2 * k1 * k2 * (1 - k2_ratio_term)
    b = delta_t * k3
    c = 1 + delta_t * k1 * z0 * (1 - k2_ratio_term)
    d = -y0 - k2_ratio_term * x0

    # Solve polynomial
    roots = np.roots([a, b, c, d])

    # y is the positive real root
    y = roots[np.logical_and(np.isreal(roots), roots > 0)][0]

    # Compute x and z
    z = z0 + delta_t * k3 * y**2
    x = (x0 + delta_t * k1 * y * z) / (1 + delta_t * k2)

    return np.array([x, y, z])

It is not necessary to use this, since fixed point iteration works fine.

---

References

• Hoggett, J. G. and Kellett, G. L., Kinetics of the monomer-dimer reaction of yeast hexokinase PI, Biochem. J., 287, 567–572, 1992. (link)

• Robertson, H. H., The soluton of a set of reaction rate equations, In Numerical Analysis: An Introduction, Walsh, J., Ed., Academic Press, 178–182, 1967.

Computing environment

%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,jupyterlab

Python implementation: CPython
Python version       : 3.10.10
IPython version      : 8.10.0

numpy     : 1.23.5
scipy     : 1.10.0
bokeh     : 3.1.0
jupyterlab: 3.5.3