Using the PMR2 generated C code for CellML models

In order to run the C code, you need to use an ODE (ordinary differential equation) IVP (initial value problem) solver.  The Sundials CVODE is a good option:

An ODE solver solves a problem that is of the form y(t)' = f(t,y(t)) with y(0) = y0.


In the generated code, f(t, y(t)) is "computeRates"


Setting y to y0 is done by calling "initConsts" (a slightly misleading name, since it also initialises state variables).

 In broad terms, an ODE solver works as follows.


You call an ODE solver's setup function to specify what function it must use as f(t,y(t)), and what it must use as initial values for y.

You then usually repeatedly call a solver time stepping function that will in turn repeatedly call f(t,y(t)). It will use the calculated rates in order to calculate the values for y at each time step.  (Note: some ODE solvers will compute a batch of time steps, rather than requiring being called for each time steps, e.g. MATLAB's ode15s and related MATLAB ODE solvers.)

The variables are as follows:


t is VOI (variable of integration, but for succinctness, this will be referred to below as "time"; but note that CellML does not have the restriction that VOI is always "time".)


y is "STATES" (an array)

y' is "RATES" (also an array)

"CONSTANTS" (also an array) are just that, constants.  It may seem that they could have been generated as literals, but CellML allows constants to be specified by expressions involving other constants, hence this approach was chosen.


As for ALGEBRAIC (also an array), it requires a bit more explanation:


Firstly, although we are usually treat the model as an ODE, most models currently in the repository are strictly speaking index 1 DAEs (differential algebraic equations).  In simple terms, this is because some of the time-varying state variables do not have their derivatives taken, and we call these "algebraic" variables.  They are also often sometimes called intermediate variables.  The reason that we can treat most of the current repository's index 1 DAEs as if they ODEs is that in these cases, we can essentially keep the ODE solver oblivious to the existence of the "algebraic" variables.  However, there is a catch: if we do this, the values that are left in the ALGEBRAIC array after the solver returns values for STATES for each time step do not usually correspond to the values in the STATES array. In other words, where an ALGEBRAIC's value depends on one or more state variables, the value in ALGEBRAICS will not match the value that would be computed using the values in STATES. This is because the solver usually uses the values calculated from several calls to f(t, y(t)) to compute the value for y(t + timestep), and does not necessarily call f(t, y(t)) again after it has calculated y(t + timestep).

Secondly, not all of the ALGEBRAIC variables are required in order to evaluate f(t, y(t)), and so to save time during integration, these are omitted from "computeRates".  Since an ODE solver will usually call "computeRates" multiple times for each time step, it would be a waste to calculate unused ALGEBRAICS inside of f(t, y(t)).


For the above two reasons, we provide "computeVariables".  It can be used after each time step to calculate correct values of the algebraic that were needed within f(t,y(t)), and simultaneously to calculate the rest of the algebraics.


If you are not interested in the values of the algebraic, you don't need to call "computeVariables"


One more thing: there are a small number of models in the CellML model repository at the moment that are index 1 DAE's where the trick of treating them as an ODE does not really work.  For these cases, it is better to use a full-blown DAE solver, and here we have had the most success with Sundials IDA. It is possible to generate IDA style code for these using the CellML API, but we have not yet made this feature available on the web interface to the repository. 


If you are dealing with a model that does fall into this category, it is currently best to obtain the CellML-API itself, which can be used to generate the required code, and will even run the simulations for you.  The easiest way of doing this is via the OpenCell software which uses the CellML-API, but it is also fairly straight forward to also do this directly from the CellML-API.