CellMLToolkit.jl
CellMLToolkit.jl is a Julia library that connects CellML models to SciML, the Scientific Julia ecosystem. CellMLToolkit.jl acts as a bridge between CellML and ModelingToolkit.jl. It imports a CellML model (in XML) and emits a ModelingToolkit.jl intermediate representation (IR), which can then enter the SciML ecosystem.
CellML
CellML is an XML-based open-standard for the exchange of mathematical models. CellML originally started in 1998 by the Auckland Bioengineering Institute at the University of Auckland and affiliated research groups. Since then, its repository has grown to more than a thousand models. While CellML is not domain-specific, its focus has been on biomedical models. Currently, the active categories in the repository are Calcium Dynamics, Cardiovascular Circulation, Cell Cycle, Cell Migration, Circadian Rhythms, Electrophysiology, Endocrine, Excitation-Contraction Coupling, Gene Regulation, Hepatology, Immunology, Ion Transport, Mechanical Constitutive Laws, Metabolism, Myofilament Mechanics, Neurobiology, pH Regulation, PKPD, Protein Modules, Signal Transduction, and Synthetic Biology. There are many software tools to import, process and run CellML models; however, these tools are not Julia-specific.
SciML
SciML is a collection of Julia libraries for open source scientific computing and machine learning. The centerpiece of SciML is DifferentialEquations.jl, which provides a rich set of ordinary differential equations (ODE) solvers. One major peripheral component of SciML is ModelingToolkit.jl. It is a modeling framework for high-performance symbolic-numeric computation in scientific computing and scientific machine learning. The core of ModelingToolkit.jl is an IR language to code the scientific problems of interest in a high level. Automatic code generation and differentiation allow for the generation of a usable model for the other components of SciML, such as DifferentialEquations.jl.
Install
To install, run
Pkg.add("CellMLToolkit")
Simple Example
using CellMLToolkit, DifferentialEquations, Plots
ml = CellModel("models/lorenz.cellml.xml")
tspan = (0, 100.0)
prob = ODEProblem(ml, tspan)
sol = solve(prob, TRBDF2(), dtmax=0.01)
X = map(x -> x[1], sol.u)
Z = map(x -> x[3], sol.u)
plot(X, Z)
Tutorial
The models directory contains 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 states 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,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)
).
Source Generation
In many applications, it is easier if we have a standalone Julia file defining our model instead of regenerating the model on the fly each time. CellMLToolkit.generate_files
compiles a CellML
model into a Julia file that can be used independent of CellMLToolkit.jl and ModelingToolkit.jl. For example, we can convert the Lorenz model by running the following code:
using CellMLToolkit
ml = CellModel("models/lorenz.cellml.xml"; dependency=false)
CellMLToolkit.generate_functions("lorenz.jl", ml)
Note, dependency=false
is needed here. This function generates a Julia source code named lorenz.jl
that contains two functions f!
(derivatives) and J!
(the Jacobian). Passing false
as the third argument of generate_functions
suppresses Jacobian generation. Also, the initial condition (u0
) and parameters (p
) are exported.
H(x) = (x >= zero(x) ? one(x) : zero(x))
# initial conditions
u0 = [1.0, 1.0, 1.0]
# parameters
p = [10.0, 2.66667, 28.0]
function f!(duₚ, uₚ, pₚ, tₚ)
t = tₚ
# state variables:
x = uₚ[1]
z = uₚ[2]
y = uₚ[3]
# parameters:
sigma = pₚ[1]
beta = pₚ[2]
rho = pₚ[3]
# algebraic equations:
# system of ODEs:
∂x = sigma * (y - x)
∂y = x * (rho - z) - y
∂z = x * y - beta * z
# state variables:
duₚ[1] = ∂x
duₚ[2] = ∂z
duₚ[3] = ∂y
nothing
end
function J!(J, uₚ, pₚ, tₚ)
t = tₚ
# state variables:
x = uₚ[1]
z = uₚ[2]
y = uₚ[3]
# parameters:
sigma = pₚ[1]
beta = pₚ[2]
rho = pₚ[3]
# Jacobian:
J[1,1] = -sigma
J[1,2] = sigma
J[2,1] = -z + rho
J[2,2] = -1
J[2,3] = -x
J[3,1] = y
J[3,2] = x
J[3,3] = -beta
nothing
end
We can run this file, independent of CellMLToolkit, as
using DifferentialEquations, Plots
include("lorenz.jl")
prob = ODEProblem(f!, u0, tspan, p)
sol = solve(prob, Tsit5(), dtmax=0.01)
sol = Array(sol)
plot(sol[1,:], sol[3,:])
The main benefit of this method is the ability to edit and modify the source code instead of the XML tree describing a CellML model.