Tutorial
The models directory contains a few CellML model examples. Let's start with a simple one, the famous Lorenz equations!
using CellMLToolkit
ml = CellModel("models/lorenz.cellml.xml")
tspan = (0, 100.0)
prob = ODEProblem(ml, tspan)
Now, ml
points to a CellModel
struct that contains the details of the model, and prob
is an ODEProblem
ready for integration. We can solve and visualize prob
as
using DifferentialEquations, Plots
sol = solve(prob, TRBDF2(), dtmax=0.01)
X = map(x -> x[1], sol.u)
Z = map(x -> x[3], sol.u)
plot(X, Z)
As expected,
Let's look at more complicated examples. The next one is the ten Tusscher-Noble-Noble-Panfilov human left ventricular action potential model. This is a mid-range electrophysiology model with 17 state variables and relatively good numerical stability.
ml = CellModel("models/tentusscher_noble_noble_panfilov_2004_a.cellml.xml")
tspan = (0, 5000.0)
prob = ODEProblem(ml, tspan)
sol = solve(prob, TRBDF2(), dtmax=1.0)
V = map(x -> x[1], sol.u)
plot(sol.t, V)
We can also enhance the model by asking ModelingToolkit.jl to generate a Jacobian by passing jac=true
to the ODEProblem
constructor.
prob = ODEProblem(ml, tspan; jac=true)
The rest remains the same. For the last example, we chose a complex model to stress the ODE solvers: the O'Hara-Rudy left ventricular model. This model has 49 state variables, is very stiff, and is prone to oscillation. The best solver for this model is CVODE_BDF
from the Sundial suite.
ml = CellModel("models/ohara_rudy_cipa_v1_2017.cellml.xml")
tspan = (0, 5000.0)
prob = ODEProblem(ml, tspan);
sol = solve(prob, CVODE_BDF(), dtmax=0.5)
V = map(x -> x[1], sol.u)
plot(sol.t, V)
Changing Parameters
Up to this point, we have run the model exactly as provided by CellML. In practice, we need to be able to modify the model parameters (either the initial conditions or the proper parameters). CellMLToolkit has multiple utility functions that help us interrogate and modify the model parameters.
There are three list
functions: list_states
, list_params
, and list_initial_conditions
. list_states
returns a list of the state variables, i.e., the variables present on the left side of an ODE. list_params
and list_initial_conditions
return arrays of (variable, value) pairs, providing the model parameters and the state variables initial conditions, respectively (corresponding to p
and u0
in DifferentialEquations.jl nomenclature).
Here, we are interested in list_params
. Let's go back to the ten Tusscher-Noble-Noble-Panfilov model and list its params:
ml = CellModel("models/tentusscher_noble_noble_panfilov_2004_a.cellml.xml")
p = list_params(ml)
display(p)
We get a list of the 45 parameters:
45-element Array{Pair{Operation,Float64},1}:
stim_start => 10.0
g_pK => 0.0146
g_bna => 0.00029
K_mNa => 40.0
b_rel => 0.25
g_Ks => 0.062
K_pCa => 0.0005
g_Kr => 0.096
Na_o => 140.0
K_up => 0.00025
⋮
To modify a parameter, we use update_list!
function. For example, the following code changes the stimulation period (stim_period
) from its default of 1000 ms to 400 ms
update_list!(p, "stim_period", 400.0)
We need to pass the new p
to ODEProblem
constructor as a keyword parameter. The rest of the code remains the same.
tspan = (0, 5000.0)
prob = ODEProblem(ml, tspan; p = p)
sol = solve(prob, TRBDF2(), dtmax = 1.0)
V = map(x -> x[1], sol.u)
plot(sol.t, V)
ODEProblem
also accepts a u0
parameter to change the initial conditions (remember u0 = list_initial_conditions(ml)
).