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: ˙y(t)=f(y(t),u(t),t,p)y(t0)=y0
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 t0 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
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.
Step sizes:
CV_ONE_STEP
approach where CVODE just makes a single step of whatever length it determines is best.Models of single-cell electrophysiology have two properties that make CVODE a good choice:
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.
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.
Myokit can use CVODES to perform forward sensitivity analysis. In this method, we calculate both the solution to our ODEs y(t)=∫tt0f(y,u,τ,p)dτ, and the partial derivatives ∂∂qy(t) of the solution with respect to some parameter or initial condition q.
In other words, given a system ˙y(t)=f(y(t),u(t),t,p)y(t0)=y0(q)p=p(q)
we want to evaluate both y(t)=∫tt0f(y(τ),u(τ),τ,p)
where y and si are vector quantities (with one entry per state variable).
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 y0 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:
˙y(t)=f(y,t,q)y(t0)=y0(q)We can now define the sensitivity si(t) as a vector quantity containing the partial derivatives of the state vector y(t) with respect to qi:
si(t)=∂∂qiy(t)In forward sensitivity analysis, we write an expression for the time derivative of the sensitivity, ˙si, and integrate this along with f to obtain si.
We find ˙si as
˙si(t)=ddt∂∂qiy(t,qi)=∂∂qi[ddty(t,qi)]Here, the expression [ddty(t,qi)] is a function of t and qi, which we can evaluate using f.
However, f is a function of t, qi, and y(qi), f=f(y(t,qi),t,qi), so that we can't simply substitute it in. Instead, we need to apply the multivariate chain rule to obtain:
˙si(t)=∂f∂y∂y∂qi+∂f∂t∂t∂qi+∂f∂qi∂qi∂qi=∂f∂y∂y∂qi+∂f∂qiwhere we recognise the Jacobian matrix ∂f∂y making an appearance, as well as the quantity ∂y∂qi=si.
˙si(t)=∂f∂ysi(t)+∂f∂qiFinally, to get the sensitivity of the j-th state in the vector, we select the j-th entry from ˙si(t) as ˙sij(t). For the first term, this amounts to selecting the dot product of row vector ∂fj∂y and column vector si(t) from the matrix multipilcation, so that we find:
˙sij(t)=∂fj∂qi+n∑k=1∂fj∂yksikLike with the Jacobian, CVODES can use finite difference approximation to estimate ∂f∂qi, and since it already knows the Jacobian and the previous sensitivities, it can calculate
˙si(t)=∂f∂ysi(t)+∂f∂qiand 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.
TODO
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]/[pi]. With this, the absolute tolerance for the j-th component of the sensitivity vector si is set to atolj/|ˉpi|, where atolj are the absolute tolerances for the state variables and ˉ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.
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 G(q)=∫tt0g(y,τ,q)dτ
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.
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.