Technical note: CVODE(S) simulation¶

Myokit's Simulation class uses the CVODE routine from the Sundials package to integrate systems of ODEs.

In these technical notes, we write the system as:

$$\begin{align} \dot{y}\left(t\right) &= f\left(y\left(t\right),u\left(t\right),t,p\right) \\ y(t_0) &= y_0 \end{align}$$

where $y$ is the state, $u$ are external inputs, and $p$ are parameters. Note that $u$ in this case consists only of the dimensionless "pacing" signal (see the technical notes on pacing) - there are no diffusion currents involved.

Using CVODE, we can estimate $y(t)$ by numerically integrating $f$ from $t_0$ to $t$.

In addition, Sundials provides a CVODES method that extends CVODE with the capability to perform forward sensitivity analysis, i.e. to calculate the partial derivatives of the states with respect to parameters or initial conditions

Documentation for CVODE and CVODES can be found on the Sundials website

CVODE and CVODES¶

Starting from version 1.33.0, Myokit uses the CVODES C api for simulations with and without sensitivities. However, as CVODES-without-sensitivities reduces to CVODE, we use both names throughout the documentation.

Implementation choices¶

• Myokit uses CVODE's Backwards Different Formulas (BDF) integrator for stiff methods.
• Because this is an implicit method, at each step a nonlinear system must be solved, for which the "Newton iteration" method (for dense systems) is chosen.
• During integration, CVODE estimates the Jacobian $\frac{\partial f}{\partial y}$ using "difference quotients" (see the CVODE/CVODES manual section 2.1).

Step sizes:

• Instead of asking CVODES to integrate from some $t_\text{min}$ to $t_\text{max}$, Myokit uses the CV_ONE_STEP approach where CVODE just makes a single step of whatever length it determines is best.
• If this step passes over discontinuity (as indicated by the pacing protocol) at some time $t$, the state at $t$ is obtained (from CVODE, which uses a higher-order interpolation), and the simulation is reset to that point.
• If this step passes over one or multiple logging points, the state at those points is estimated (using interpolation again) without resetting the simulator.

Why CVODE(s)?¶

Models of single-cell electrophysiology have two properties that make CVODE a good choice:

• They are stiff, containing processes that occur at very different time scales. E.g. INa activates in 1-2ms, while concentrations in the same model drift very slowly. Stiff problems are difficult to integrate, and usually require an implicit ODE solver. Similarly, because the fast processes are only active for brief periods during the simulation, an adaptive step-size method is highly desirable.
• Models of single cell electrophysiology commonly have an expensive-to-evaluate right-hand side (RHS). This means that schemes requiring several RHS evaluations per step may not be competitive.

CVODE uses a variable-order (max 5) BDF integrator, which has good stability properties, so that it can take large steps even for stiff problems. You can try this out: Run a simulation of a modern cell model in CVODES, and store every step. Next, use a forward-Euler method with the same steps. When I tried this the Euler method produced highly inaccurate results. CVODE uses a predictor-corrector method to estimate its error at each step, and choose an optimal step size.

A second advantage of the BDF methods is that they are multistep methods. This means that, instead of requiring several iterations per time step, the method uses the derivatives calculated during previous iterations to build up a higher-order estimate. In Clerx & Collins 2014 we found that CVODE required only 1 or 2 RHS evaluations for most steps, although this occasionally rose dramatically, leading to an average number of evaluations per step of 1.6 to 6 for the models tested.

Evaluating the Jacobian¶

Instead of letting CVODE(s) estimate the Jacobian, it is also possible to write a "user-provided" function that calculates it.

Options for doing this automatically include automatic differentiation (this is included in Myokit but possibly won't be for long) and symbolic differentiation. Cooper et al. tried out a symbolically generated Jacobian in Chaste (which also uses CVODE) and found moderate speed-ups for most models. They also indicated a downside: the equations for the jacobian can contain singularities the same way the normal RHS does, which can lead to simulation problems.

Forward sensitivity analysis¶

Myokit can use CVODES to perform forward sensitivity analysis. In this method, we calculate both the solution to our ODEs $y(t) = \int_{t_0}^t f(y, u, \tau, p) d\tau$, and the partial derivatives $\frac{\partial}{\partial q}y(t)$ of the solution with respect to some parameter or initial condition $q$.

In other words, given a system

$$\begin{align} \dot{y}(t) &= f\left(y(t), u(t), t, p\right) \\ y(t_0) &= y_0(q) \\ p &= p(q) \end{align}$$

we want to evaluate both

$$\begin{equation} y(t) = \int_{t_0}^{t} f\left(y(\tau), u(\tau), \tau, p\right) \end{equation}$$ and $$\begin{equation} s_i(t) = \frac{\partial}{\partial q_i} y(t) \end{equation}$$

where $y$ and $s_i$ are vector quantities (with one entry per state variable).

Notation¶

To simplify things, and stay closer to the notation used in the CVODES manual, we now introduce two slight changes in notation.

First, we stop explicitly mentioning the forcing term $u(t)$, so that we can write $f(y, t, q)$. Since we want to look at partial derivatives, this is an acceptable shorthand as long as $u$ does not depend on $q$ or $y(q)$. The pacing system is already defined to be independent of $y$, so the second condition is automatically met (see also the technical notes on pacing). Having $u$ be independent of $q$ means that we cannot calculate sensitivities w.r.t. the pacing protocol, which may be a downside for protocol-design projects, but is very easy to implement. Having $u$ independent of $y$ means that $u$ can not include diffusion currents in this analysis. As we are concerned with single-cell simulations in this case, this condition is also easily met.

We will deviate from the CVODES manual in our use of $p$: In the CVODES formulation this is a "parameter vector" that both $f$ and $y_0$ depend on. In Myokit models, the initial state cannot depend on any model parameters, and so we use $q$ to indicate a vector consisting of selected model parameters and initial conditions. Next, we define $p$ to only include parameters independent of $q$. In other words $p$ contains the parameters not used in sensitivity analysis, while $q$ contains the parameters and initial conditions used in sensitivity analysis.

Just like $u(t)$, we can then omit $p$ from our notation, so that the system is written as:

$$\begin{align} \dot{y}(t) &= f(y, t, q) \\ y(t_0) &= y_0(q) \end{align}$$

Sensitivities¶

We can now define the sensitivity $s_i(t)$ as a vector quantity containing the partial derivatives of the state vector $y(t)$ with respect to $q_i$:

$$\begin{equation} s_i(t) = \frac{\partial}{\partial q_i} y(t) \end{equation}$$

In forward sensitivity analysis, we write an expression for the time derivative of the sensitivity, $\dot{s}_i$, and integrate this along with $f$ to obtain $s_i$.

We find $\dot{s}_i$ as

$$\begin{equation} \dot{s}_i(t) = \frac{d}{d t} \frac{\partial}{\partial q_i} y(t, q_i) = \frac{\partial}{\partial q_i} \left[ \frac{d}{d t} y(t, q_i) \right] \end{equation}$$

Here, the expression $\left[ \frac{d}{d t} y(t, q_i) \right]$ is a function of $t$ and $q_i$, which we can evaluate using $f$.

However, $f$ is a function of $t$, $q_i$, and $y(q_i)$, $f = f\left(y(t, q_i), t, q_i\right)$, so that we can't simply substitute it in. Instead, we need to apply the multivariate chain rule to obtain:

$$\begin{align} \dot{s}_i(t) &= \frac{\partial f}{\partial y} \frac{\partial y}{\partial q_i} + \frac{\partial f}{\partial t} \frac{\partial t}{\partial q_i} + \frac{\partial f}{\partial q_i} \frac{\partial q_i}{\partial q_i} \\ &= \frac{\partial f}{\partial y} \frac{\partial y}{\partial q_i} + \frac{\partial f}{\partial q_i} \end{align}$$

where we recognise the Jacobian matrix $\frac{\partial f}{\partial y}$ making an appearance, as well as the quantity $\frac{\partial y}{\partial q_i} = s_i$.

$$\begin{align} \dot{s}_i(t) = \frac{\partial f}{\partial y} s_i(t) + \frac{\partial f}{\partial q_i} \end{align}$$

Finally, to get the sensitivity of the j-th state in the vector, we select the $j$-th entry from $\dot{s}_i(t)$ as $\dot{s}_{ij}(t)$. For the first term, this amounts to selecting the dot product of row vector $\frac{\partial f_j}{\partial y}$ and column vector $s_i(t)$ from the matrix multipilcation, so that we find:

$$\begin{align} \dot{s}_{ij}(t) = \frac{\partial f_j}{\partial q_i} + \sum_{k=1}^{n} \frac{\partial f_j}{\partial y_k} s_{ik} \end{align}$$

Implementation (approximate)¶

Like with the Jacobian, CVODES can use finite difference approximation to estimate $\frac{\partial f}{\partial q_i}$, and since it already knows the Jacobian and the previous sensitivities, it can calculate

$$\begin{equation} \dot{s}_i(t) = \frac{\partial f}{\partial y} s_i(t) + \frac{\partial f}{\partial q_i} \end{equation}$$

and integrate this along with the other ODEs in the system.

What's left to do is to calculate the sensitivities of intermediate variables with respect to $q$. This is hanlded by an extra step that uses symbolic differentiation to calculate the intermediate sensitivities based on the current state-sensitivities, which is called before any logging operation.

Implementation (symbolic)¶

TODO

Error control¶

For error control, CVODES uses absolute and relative tolerances:

The relative tolerance for sensitivity variables is set to be the same as for the state variables. The selection of absolute tolerances for the sensitivity variables is based on the observation that the sensitivity vector s i will have units of $[y]/[p_i]$. With this, the absolute tolerance for the $j$-th component of the sensitivity vector $s_i$ is set to $\text{atol}_j / |\bar{p}_i|$, where $\text{atol}_j$ are the absolute tolerances for the state variables and $\bar{p}$ is a vector of scaling factors that are dimensionally consistent with the model parameters $p$ and give an indication of their order of magnitude.

Note that their $p$ is $q$ in our notation. At the moment, Myokit sets these scaling factors based on the values of $q$.

The "Simultaneous Corrector" mode of error correction was chosen (CVODES 2.6.1) as it seemed the most appropriate. No experiments to verify this have been performed.

Adjoint sensitivity analysis¶

An alternative method for sensitivity analysis that is supported in CVODES is adjoint sensitivity analysis.

In adjoint sensitivity analysis the user has some function g(y, t, q) that gets integrated over several time points, for which a gradient must be found.

So given

$$\begin{equation} G(q) = \int_{t_0}^{t} g(y, \tau, q) d\tau \end{equation}$$ we want to find $\frac{dG}{dp}$.

This form might be useful for parameter estimation, in which case $G$ would be a score function or a likelihood.

However, this would require the user to be able to supply an equation (in some format myokit understands) for $g$, which doesn't fit with the current goal of moving all parameter estimation methods in Myokit to PINTS.

CModel class¶

In Myokit version 1.33.0 and after, an attempt has been made to separate the CVODE(s) simulator code from the model code wherever possible. This is described in detail in the cmodel.h header template, which implements a C model interface for Myokit models.