Background to the VPH-Physiome project
To be of benefit to applications in healthcare, organ and whole organism
physiology needs to be understood at both a systems level and in terms
of subcellular function and tissue properties. Understanding a
re-entrant arrhythmia in the heart, for example, depends on knowledge of
not only numerous cellular ionic current mechanisms and signal
transduction pathways, but also larger scale myocardial tissue structure
and the spatial variation in protein expression. As reductionist
biomedical science succeeds in elucidating ever more detail at the
molecular level, it is increasingly difficult for physiologists to
relate integrated whole organ function to underlying biophysically
detailed mechanisms that exploit this molecular knowledge. Multi-scale
computational modelling is used by engineers and physicists to design
and analyse mechanical, electrical and chemical engineering systems.
Similar approaches could benefit the understanding of physiological
systems. To address these challenges and to take advantage of
bioengineering approaches to modelling anatomy and physiology, the
International Union of Physiological Sciences (IUPS) formed the Physiome
Project in 1997 as an international collaboration to provide a
computational framework for understanding human physiology .
Primary Goals
One of the primary goals of the Physiome Project [PJ04] has been to promote
the development of standards for the exchange of information between
models. The first of these standards, dealing with time varying but
spatially lumped processes, is CellML [VarYY]. The second (dealing with
spatially and time varying processes) is FieldML [CPJ09][P13] . A further
goal of the Physiome Project has been the development of open source
tools for creating and visualizing standards-based models and running
model simulations. OpenCOR is the latest in a series of software
projects aimed at providing a modelling environment for CellML models.
Similar tools exist for FieldML models.
Following the publication of the STEP (Strategy for a European
Physiome) Roadmap in 2006, the European Commission in 2007 initiated
the Virtual Physiological Human (VPH) project [ea13]. A related US
initiative by the Interagency Modeling and Analysis Group (IMAG) began
in 2003 . These projects and similar initiatives are now coordinated
and are collectively referred to here as the ‘VPH-Physiome’
project . The VPH-Institute was formed in 2012 as a virtual
organisation to providing strategic leadership, initially in Europe but
now globally, for the VPH-Physiome Project.
Footnotes
Create and run a simple CellML model: editing and simulation
In this example we create a simple CellML model and run it. The model is
the Van der Pol oscillator defined by the second order equation
\[\frac{d^{2}x}{dt^{2}} - \mu\left( 1 - x^{2} \right)\frac{\text{dx}}{\text{dt}} + x = 0\]
with initial conditions
\(x = - 2;\ \frac{\text{dx}}{\text{dt}} = 0\). The parameter
\(\mu\) controls the magnitude of the damping term. To create a
CellML model we convert this to two first order equations by
defining the velocity \(\frac{\text{dx}}{\text{dt}}\) as a new
variable \(y\):
()\[\frac{\text{dx}}{\text{dt}} =\ y\]
()\[\frac{\text{dy}}{\text{dt}} =\ \mu\left( 1 - x^{2} \right)y - x\]
The initial conditions are now \(x = - 2;y = 0\).
With the central pane in Editing mode (e.g. CellML Text view), create a new CellML file:
and then type in the
following lines of code after deleting the three lines that indicate
where the code should go:
def model van_der_pol_model as
def comp main as
var t: dimensionless {init: 0};
var x: dimensionless {init: -2};
var y: dimensionless {init: 0};
var mu: dimensionless {init: 1};
// These are the ODEs
ode(x,t)=y;
ode(y,t)=mu*(1{dimensionless}-sqr(x))*y-x;
enddef;
enddef;
Things to note are:
- the closing semicolon at the end of each line (apart from the first two def statements that are opening a CellML construct);
- the need to indicate dimensions for each variable and constant (all dimensionless in this example – but more on dimensions later);
- the use of ode(x,t) to indicate a first order ODE in x and t
- the use of the squaring function sqr(x) for \(x^{2}\), and
- the use of ‘//’ to indicate a comment.
A partial list of mathematical functions available for OpenCOR is:
| \(x^{2}\) |
sqr(x) |
\(\sqrt{x}\) |
sqrt(x) |
\(\ln x\) |
ln(x) |
\(\log_{10}x\) |
log(x) |
\(e^{x}\) |
exp(x) |
\(x^{a}\) |
pow(x,a) |
| \(\sin x\) |
sin(x) |
\(\cos x\) |
cos(x) |
\(\tan x\) |
tan(x) |
\(\csc x\) |
csc(x) |
\(\sec x\) |
sec(x) |
\(\cot x\) |
cot(x) |
| \(\sin^{-1}x\) |
asin(x) |
\(\cos^{-1} x\) |
acos(x) |
\(\tan^{-1} x\) |
atan(x) |
\(\csc^{-1} x\) |
acsc(x) |
\(\sec^{-1} x\) |
asec(x) |
\(\cot^{-1}x\) |
acot(x) |
| \(\sinh x\) |
sinh(x) |
\(\cosh x\) |
cosh(x) |
\(\tanh x\) |
tanh(x) |
\(\operatorname{csch} x\) |
csch(x) |
\(\operatorname{sech} x\) |
sech(x) |
\(\coth x\) |
coth(x) |
| \(\sinh^{-1} x\) |
asinh(x) |
\(\cosh^{-1} x\) |
acosh(x) |
\(\tanh^{-1} x\) |
atanh(x) |
\(\operatorname{csch}^{-1} x\) |
acsch(x) |
\(\operatorname{sech}^{-1}x\) |
asech(x) |
\(\coth^{-1} x\) |
acoth(x) |
Table 1. Partial list of mathematical functions available for coding in
OpenCOR.
Positioning the cursor over either of the ODEs renders the maths in
standard form above the code as shown in Fig. 1.
Note that CellML is a declarative language (unlike say C, Fortran
or Matlab, which are procedural languages) and therefore the order of
statements does not affect the solution. For example, the order of the
ODEs could equally well be
ode(y,t)=mu*(1{dimensionless}-sqr(x))*y-x;
ode(x,t)=y;
The significance of this will become apparent later when we import
several CellML models to create a composite model.
Now save the code to a local folder using Save under the File menu ()
(or ‘CTRL-S’) and choosing .cellml as the file format. With the
CellML model saved various views, accessed via the tabs on the right
hand edge of the window, become available. One is the CellML Text view
(the view used to enter the code above); another is the Raw CellML
view that displays the way the model is stored and is intentionally
verbose to ensure that the meaning is always unambiguous (note that
positioning the cursor over part of the code shows the maths in this
view also); and another is the Raw view. Notice that ‘CTRL-T’ in the
Raw CellML view performs validation tests on the CellML model. The
CellML Text view provides a much more convenient format for entering
and editing the CellML model.
With the equations and initial conditions defined, we are ready to run
the model. To do this, click on the Simulation tab on the left hand
edge of the window. You will see three main areas - at the left hand
side of the window are the Simulation, Solvers, Graphs and
Parameters panels, which are explained below. At the right hand side
is the graphical output window, and running along the bottom of the
window is a status area, where status messages are displayed.
Solvers Panel
This area is used to configure the solver that will run the simulation.
- Name - this is used to set the solver algorithm. It will be set by
default to be the most appropriate solver for the equations you are
solving. OpenCOR allows you to change this to another solver
appropriate to the type of equations you are solving if you choose
to. For example, CVODE for ODE (ordinary differential equation)
problems, IDA for DAE (differential algebraic equation) problems,
KINSOL for NLA (non-linear algebraic) problems.
- Other parameters for the chosen solver – e.g. Maximum step,
Maximum number of steps, and Tolerance settings for CVODE and
IDA. For more information on the solver parameters, please refer to
the documentation for the particular solver.
Note: these can all be left at their default values for our simple demo
problem.
Graphs Panel
This shows what parameters are being plotted once these have been
defined in the Parameters panel. These can be selected/deselected by
clicking in the box next to a parameter.
Parameters Panel
This panel lists all the model parameters, and allows you to select one
or more to plot against the variable of integration or another parameter
in the graphical output windows. OpenCOR supports graphing of any
parameter against any other. All variables from the model are listed
here, arranged by the components in which they appear, and in
alphabetical order. Parameters are displayed with their variable name,
their value, and their units. The icons alongside them have the
following meanings:
 Editable constant |
 Editable state variable |
 Computed constant |
 Rate variable |
 Variable of integration |
 Algebraic quantity |
Right clicking on a parameter provides the options for displaying that
parameter in the currently selected graphical output window. With the
cursor highlighting the top graphical output window (a blue line appears
next to it), select x then Plot Against Variable of Integration – in
this case t - in order to plot x(t). Now move the cursor to the
second graphical output window and select y then t to plot y(t).
Finally select the bottom graphical output window, select y and select
Plot Against then Main then x to plot y(x).
Now click on the Run control. You will see a progress bar running
along the bottom of the status window. Status messages about the
successful simulation, including the time taken, are displayed in the
bottom panel. This can be hidden by dragging down on the bar just above
the panel. Fig. 2 shows the results. Use the interval delay wheel to
slow down the plotting if you want to watch the solution evolve. You can
also pause the simulation at any time by clicking on the Run control
and if you change a parameter during the pause, the simulation will
continue (when you click the Run control button again) with the new
parameter.
Note that the values shown for the various parameters are the values
they have at the end of the solution run. To restore these to their
initial values, use the Reset parameters (
) button. To clear
the graphical output traces, click on the Clear simulation data
(
) button.
The top two graphical output panels are showing the time-dependent
solution of the x and y variables. The bottom panel shows how y
varies as a function of x. This is called the solution in state space
and it is often useful to analyse the state space solution to capture
the key characteristics of the equations being solved.
To obtain numerical values for all variables (i.e. x(t) and y(t)),
click on the CSV file button (
). You will be asked to enter a
filename and type (use .csv). Opening this file (e.g. with Microsoft
Excel) provides access to the numerical values. Other output types (e.g.
BiosignalML) will be available in future versions of OpenCOR.
You can move the graphical output traces around with ‘left click and
drag’ and you can change the horizontal or vertical scale with ‘right
click and drag’. Holding the SHIFT key down while clicking on a
graphical output panel allows you to interrogate the solution at any
point. Right clicking on a panel provides zoom facilities.
Note
The simulation described above can also be loaded and run directly in OpenCOR using this link.
The various plugins used by OpenCOR can be viewed under the Tools menu.
A French language version of OpenCOR is also available under the Tools
menu. An option under the File menu allows a file to be locked (also
‘CTRL-L’). To indicate that the file is locked, the background colour
switches to pink in the CellML Text and Raw CellML views and a
lock symbol appears on the filename tab. Note that OpenCOR text is case
sensitive.
Footnotes
Open an existing CellML file from a local directory or the Physiome Model Repository
Go to the File menu and select Open… (). Browse to the folder that
contains your existing models and select one. Note that this brings up a
new tabbed window and you can have any number of CellML models open at
the same time in order to quickly move between them. A model can be
removed from this list by clicking on 
next to the CellML model
name.
You can also access models from the left hand panel in Fig. 1.1(a). If
this panel is not currently visible, use ‘CTRL-spacebar’ to make it
reappear. Models can then be accessed from any one of the three
subdivisions of this panel – File Browser, Physiome Model Repository
or File Organiser. For a file under File Browser or File
Organiser, either double-click it or ‘drag&drop’ it over the central
workspace to open that model. Clicking on a model in the Physiome Model
Repository (PMR) (e.g. Chen, Popel, 2007) opens a new browser window
with that model (PMR is covered in more detail in Section 13). You can
either load this model directly into OpenCOR or create an identical copy
(clone) of the model in your local directory. Note that PMR contains
workspaces and exposures. Workspaces are online environments for the
collaborative development of models (e.g. by geographically dispersed
groups) and can have password protected access. Exposures are workspaces
that are exposed for public view and mostly contain models from
peer-reviewed journal publications. There are about 600 exposures based
on journal papers and covering many areas of cell processes and other
ODE/algebraic models, but these are currently being supplemented with
reusable protein-based models – see discussion in a Section 13.
To load a model directly into OpenCOR, click on the right-most of the
two buttons in Fig. 3 - this lists the CellML models in that exposure
- and then click on the model you want. Clicking on the left hand button
copies the PMR workspace to a local directory that you specify. This is
useful if you want to use that model as a template for a new one you are
creating.
In the PMR window (Fig. 3) the buttons on the right-hand side [1] lists all the CellML files for this model. Clicking on one of those [2] uploads the model into OpenCOR. The left-hand buttons [3] copies the PMR workspace to a local directory.
A simple first order ODE
The simplest example of a first order ODE is
\[\frac{\text{dy}}{\text{dt}} = - ay + b\]
with the solution
\[y\left( t \right) = \frac{b}{a} + \left( y\left( 0 \right) - \frac{b}{a} \right).e^{- at},\]
where \(y\left( 0 \right)\) or \(y_{0}\), the value of
\(y\left( t \right)\) at \(t = 0\), is the initial condition.
The final steady state solution as \(t \rightarrow \infty\) is
\(y\left( \left. \ t \right|_{\infty} \right) = y_{\infty} = \frac{b}{a}\)
(see Figure 6). Note that \(t = \tau = \frac{1}{a}\) is called the
time constant of the exponential decay, and that
\[y\left( \tau \right) = \frac{b}{a} + \left( y\left( 0 \right) - \frac{b}{a} \right).e^{- 1}.\]
At \(t = \tau\) , \(y\left( t \right)\) has therefore fallen to
\(\frac{1}{e}\) (or about 37%) of the difference between the initial
(\(y\left( 0 \right)\)) and final steady state (
\(y\left( \infty \right)\)) values.
Choosing parameters \(a = \tau = 1;b = 2\) and
\(y\left( 0 \right) = 5\), the CellML Text for this model is
def model first_order_model as
def comp main as
var t: dimensionless {init: 0};
var y: dimensionless {init: 5};
var a: dimensionless {init: 1};
var b: dimensionless {init: 2};
ode(y,t)=-a*y+b;
enddef;
enddef;
The solution by OpenCOR is shown in Fig. 5(a) for these parameters (a
decaying exponential) and in Fig. 5(b) for parameters
\(a = 1;b = 5\) and \(y\left( 0 \right) = 2\) (an inverted
decaying exponential). Note the simulation panel with Ending
point=10, Point interval=0.1. Try putting \(a = - 1\).
These two solutions have the same exponential time constant
(\(\tau = \frac{1}{a} = 1\)) but different initial and final (steady
state) values.
The exponential decay curve shown on the left in Fig. 5 is a common
feature of many models and in the case of radioactive decay (for
example) is a statement that the rate of decay
(\(- \frac{\text{dy}}{\text{dt}}\)) is proportional to the
current amount of substance (\(y\)). This is illustrated on
the NZ$100 note (should you be lucky enough to possess one), shown in
Figure 8.
Footnotes
The Lorenz attractor
An example of a third order ODE system (i.e. three 1st order
equations) is the Lorenz equations.
This system has three equations:
\[\begin{split}\frac{\text{dx}}{\text{dt}} & = \sigma\left( y - x \right) \\
\frac{\text{dy}}{\text{dt}} & = x\left( \rho - z \right) - y \\
\frac{\text{dz}}{\text{dt}} & = xy - \beta z\end{split}\]
where \(\sigma,\ \rho\) and \(\beta\) are parameters.
The CellML Text code entered for these equations is shown in Fig. 7 with parameters
\(\sigma = 10\), \(\rho = 28\), \(\beta = 8/3\) = 2.66667
and initial conditions
\(x\left( 0 \right) = y\left( 0 \right) = z\left( 0 \right) =\)1.
Solutions for \(x\left( t \right)\), \(y\left( x \right)\) and
\(z\left( x \right)\), corresponding to the time integration
parameters shown on the LHS, are shown in Fig. 8. Note that this
system exhibits ‘chaotic dynamics’ with small changes in the initial
conditions leading to quite different solution paths.
This example illustrates the value of OpenCOR’s ability to plot
variables as they are computed. Use the Simulation Delay wheel to slow
down the plotting by a factor of about 5-10,000 - in order to follow the
solution as it spirals in ever widening trajectories around the left
hand wing of the attractor before coming close to the origin that then
sends it off to the right hand wing of the attractor.
Solutions to the Lorenz equations are organised by the 2D ‘Lorenz
manifold’. This surface has a very beautiful shape and has become an art
form - even rendered in crochet! (See Fig. 9).
Note
The simulation presented in Fig. 8 can be loaded direction into OpenCOR using this link.
Exercise for the reader
Another example of intriguing and unpredictable behaviour from a simple
deterministic ODE system is the ‘blue sky catastrophe’ model [JH02] defined
by the following equations:
\[\begin{split}\frac{\text{dx}}{\text{dt}} & = y \\
\frac{\text{dy}}{\text{dt}} & = x - x^{3} - 0.25y + A\sin t\end{split}\]
with parameter \(A = 0.2645\) and initial conditions
\(x\left( 0 \right) = 0.9\), \(y\left( 0 \right) = 0.4\). Run to
\(t = 500\) with \(\Delta t = 0.01\) and plot
\(x\left( t \right)\) and \(y\left( x \right)\) (OpenCOR link). Also try with
\(A = 0.265\) to see how sensitive the solution is to small changes
in parameter values.
Footnotes
Code generation
It is sometimes required to export CellML models to various procedural formats to make use of a given model with existing tools. OpenCOR currently uses the CellML Language Export Definition Service provided by the CellML API to achieve this (see this article for details). This service takes an XML file containing a conversion definition and uses that to export a CellML model to the defined format.
The OpenCOR distribution packages include definition files for C, Fortran 77, Python, and Matlab. These definition files are available in the formats folder of your OpenCOR installation or can be downloaded and used directly using the previous links.
The C and Fortran code generated using these definition files contain functions suitable for inclusion in DAE/ODE simulation codes. Whereas the Python and Matlab code generated are complete scripts that use standard Python or Matlab methods to actually perform an default simulation. The default simulation is probably not what is needed, so the generated code can be modified or reused to meet the specific usage requirements.
Exporting CellML to code
The steps to generate code from OpenCOR are given below.
- Load the desired CellML model into OpenCOR (both CellML 1.0 and 1.1 models can be used)
- From the OpenCOR menu, choose .
- The first file selection dialog is to provide the conversion definition file (as above).
- The second file selection dialog is to provide the file to save the generated code to.
This conversion can also be performed using OpenCOR as a command line client. In this case the command is:
$ ./OpenCOR -c CellMLTools::export myfile.cellml myformat.xml
or for a remote model:
$ ./OpenCOR -c CellMLTools::export http://mydomain.com/myfile.cellml myformat.xml
where myformat.xml can be one of the standard definition files described above.
Generated code in PMR
The Physiome Model Repository uses the same code generation service from the CellML API to generate code in the above formats for all exposures containing CellML models. These are available from the Generated Code view for CellML models. See here for an example.
Model annotation
One of the most powerful features of CellML is its ability to import
models. This means that complex models can be built up by combining
previously defined models. There is a potential problem with this
process, however, since the imported models (often developed by
completely different modellers) may represent the same biological or
biophysical entity with different expressions. The potassium channel
model in A model of the potassium channel: Introducing CellML components and connections, for example, represents the intracellular
concentration of potassium as ‘Ki’ (see the CellML Text code Potassium_ion_channel.cellml) but another model involving the intracellular potassium
concentration may use a different expression.
The solution to this dilemma is to annotate the CellML variables with
names from controlled vocabularies that have been agreed upon by the
relevant scientific community. In this case we may simply want to
annotate Ki as ‘the concentration of potassium in the cytosol’.
This expression, however, refers to three distinct entities:
concentration, potassium and cytosol. We might also want to
specify that we are referring to the cytosol of a neuron … and that the
neuron comes from a particular part of a giant squid (the experimental
animal used by Hodgkin and Huxley). Annotations can clearly get very
complicated!
What comes to our rescue here is that most scientific communities have
developed controlled vocabularies together with the relationships
between the terms of that vocabulary – called ontologies.
Furthermore relationships can always be expressed in the form
subject-predicate-object. E.g. Ki
is-the-concentration-of potassium is one relationship and
potassium in-the cytosol is another. Each object can become
the subject of another expression. We could continue, for example, with
cytosol of-the neuron, neuron of-the squid and so
on. The terms s-the-concentration-of, in-the and of-the are
the predicates and these semantically rich expressions too have to come
from controlled vocabularies. Each of these
subject-predicate-object expressions is called an RDF triple
and the World Wide Web consortium has established a framework
called the Resource Description Framework (RDF ) to support
these.
CellML models therefore contain two parts, one dealing with syntax
(the MathML definition of the models together with the structure of
components, connections, groups, units, etc) as discussed in previous
sections, and one dealing with semantics (the meanings of the
terms used in the models) discussed in this section . This latter
is also referred to as metadata – i.e. data about data.
In the CellML metadata specification the first RDF subject of a
triple is a CellML element (e.g. a variable such as ‘Ki’), the RDF
predicate is chosen from the Biomodels Biological Qualifiers
list, and the RDF object is a URI (the string of characters used to
identify the name of a resource ). Establishing these RDF links to
biological and biophysical meaning is the goal of annotation.
Note the different types of subject/object used in the RDF triples: the
concentration is a biophysical entity, potassium is a chemical
entity, the cytosol is an anatomical entity. In fact, to cover all the
terminology used in the models, CellML uses five separate ontologies:
These ontologies are available through OpenCOR’s annotation facilities
as explained below.
If we now go back to the potassium ion channel CellML model and, under
Editing, click on CellML Annotation, the various elements of the
model (Units, Components, Variables, Groups and Connections) are
displayed (see Fig. 10). If you right click on any of them a popup
menu will appear, which you can use to expand/collapse all the child
nodes, as well as remove the metadata associated with the current CellML
element or the whole CellML file. Expanding Components lists all the
components and their variables. To annotate the potassium channel
component, select it and specify a Qualifier from the list displayed:
bio:encodes, bio:isPropertyOf
bio:hasPart, bio:isVersionOf
bio:hasProperty, bio:occursIn
bio:hasVersion, bio:hasTaxon
bio:is, model:is
bio:isDescribedBy, model:isDerivedFrom
bio:isEncodedBy, model:isDescribedBy
bio:isHomologTo, model:isInstanceOf
bio:isPartOf, model:hasInstance
If you do not know which qualifier to use, click on the
button to get some information about the current qualifier
(you must be connected to the internet) and go through the list of
qualifiers until you find the one that best suits your needs. Here, we
will say that you want to use bio:isVersionOf. Fig. 11 shows the
information displayed about this qualifier.
Now you need to retrieve some possible ontological terms to describe the
potassium_channel component. For this you must enter a search term,
which in our case is ‘potassium channel’ (note that regular expressions
are supported ). This returns 24 possible ontological terms as
shown in Fig. 12. The voltage-gated potassium channel complex is the
most appropriate. Clicking on the GO identifier link shown provides more
information about this term (see Fig. 13).
Now, assuming that you are happy with your choice
of ontological term, you can associate it with the potassium_channel
component by clicking on its corresponding 
button which then displays
the qualifier, resource and ID information in the middle panel as shown
in Fig. 12. If you make a mistake, this can be removed by clicking on
the 
button.
The first level annotation of the potassium_channel component has now
been achieved. The content of the three terms in the RDF triple are
shown in Fig. 14, along with the annotation for the variables Ki and
Ko.
def comp {id_000000001} potassium_channel as
var V: millivolt {pub: in, priv: out};
var t: millisec {pub: in, priv: out};
var n: dimensionless {priv: in};
var i_K: microA_per_cm2 {pub: out};
var g_K: milliS_per_cm2 {init: 36};
var {id_000000002} Ki: mM {init: 90};
var {id_000000003} Ko: mM {init: 3};
var RTF: millivolt {init: 25};
var E_K: millivolt;
var K_conductance: milliS_per_cm2 {pub: out};
E_K = RTF*ln(Ko/Ki);
K_conductance = g_K*pow(n, 4{dimensionless});
i_K = K_conductance*(V-E_K);
enddef;
When saved (the CellML Annotation tag will appear un-grayed), the
result of these annotations is to add metadata to the CellML file. If
you switch to the CellML Text view you will see that the elements that
have been annotated appear with ID numbers, as shown above.
These point to the corresponding metadata contained in the CellML file
for this model and are displayed under the qualifier-resource-Id
headings in the annotation window when you click on the element in the
editing window.
Note that the three annotations added above are all biological
annotations. Many of the other components and variables in the CellML
potassium channel model deal with biophysical entities and these require
the use of the OPB ontology (yet to be implemented in OpenCOR). The use
of composite annotations is also being developed , such as
“Ki is-the concentration of potassium in-the
cytosol of-the neuron of-the giant-squid”, where concentration,
potassium, cytosol, neuron and giant-squid are defined by the
ontologies OPB, ChEBI, GO, FMA and a species ontology, respectively.
Footnotes
SED-ML, functional curation and Web Lab
In the same way that CellML models can be defined unambiguously, and shared easily, in a machine- readable format, there is a need to do the same thing with ‘protocols’ - i.e. to define what you have to do to replicate/simulate an experiment, and to analyse the results. An XML standard for this called SED-ML is being developed by the CellML/SBML community and preliminary support for SED-ML has been implemented in OpenCOR in order to allow precise and reproducible control over the OpenCOR simulation and graphical output (e.g., see Fig. 3.22).
The current snapshot release (2016-02-27) of OpenCOR supports exporting the Simulation view configuration to a SED-ML file, which can then be read back into OpenCOR to reproduce a given simulation experiment, illustrated in Fig. 15.
Support for SED-ML will also facilitate the curation of models according to their functional behaviour under a range of experimental scenarios.
The key idea behind functional curation is that, when mathematical and computational models are being developed, a primary goal should be the continuous comparison of those models against experimental data. When computational models are being re-used in new studies, it is similarly important to check that they behave appropriately in the new situation to which you’re applying them. To achieve this goal, a pre-requisite is to be able to replicate in-silico precisely the same protocols used in an experiment of interest. A language for describing rich ‘virtual experiment’ protocols and software for running these on compatible models is being developed in the Computational Biology Group at Oxford University.
An online system called Web Lab is also being developed that supports definition of experimental protocols for cardiac electrophysiology, and allows any CellML model to be tested under these protocols [CJ15]. This enables comparison of the behaviours of cellular models under different experimental protocols: both to characterise a model’s behaviour, and comparing hypotheses by seeing how different models react under the same protocol (Fig. 16 adapted from [CJ15]).
The Web Lab website provides tools for comparing how two different cardiac electrophysiology models behave under the same experimental protocols. Note that Web Lab demonstration for CellML models of cardiac electrophysiology is a prototype for a more general approach to defining simulation protocols for all CellML models.
Footnotes
Speed comparisons with MATLAB
Solution speed is important for complex computational models and here we
compare the performance of OpenCOR with MATLAB. Nine
representative CellML models were chosen from the PMR model repository.
For the MATLAB tests we used the MATLAB code, generated automatically
from CellML, that is available on the PMR site. These comparisons are
based on using the default solvers (listed below) available in the two
packages.
Testing environment
- MacBook Pro (Retina, Mid 2012).
- Processor: 2.6 GHz Intel Core i7.
- Memory: 16 GB 1600 MHz DDR3.
- Operating system: OS X Yosemite 10.10.3.
- Version: 0.4.1.
- Solver: CVODE with its default settings, except for its Maximum step
parameter, which is set to the model’s stimulation duration, if
needed.
- Version: R2013a.
- Solver: ode15s (i.e. a solver suitable for stiff problems and which
has low to medium order of accuracy) with both its RelTol and
AbsTol parameters set to 1e-7 and its MaxStep parameter set to
the stimulation duration, if needed.
Testing protocol
- Run a model for a given simulation duration.
- Generate simulation data every milliseconds.
- Only keep track of all the simulation data (i.e. no graphical
output).
- Run a model 7 times, discard the 2 slowest runs (to account for
unpredictable slowdowns of the testing machine) and average the
resulting computational times.
- Computational times are obtained directly from OpenCOR and MATLAB
(through a couple of calls to cputime in the case of MATLAB).
Results
*The value of membrane.stim_end was increased so as to get
action potentials for the duration of the simulation
Conclusions
For this range of tests, OpenCOR is between 38 and 218 times faster than MATLAB.
A more extensive evaluation of these results is available on GitHub.
Footnotes
Future developments
Both CellML and OpenCOR are continuing to be developed. These notes will
be updated to reflect new features of both. The next release of OpenCOR
(0.5) will include
- the SED-ML API which means that all the variables controlling the
simulation and its output can be specified in a file for that
simulation
- the BioSignalML API which will allow experimental data to be read
into OpenCOR in a standardised way
- colour plots, to better distinguish overlapping traces in the output
windows
Priorities for later releases of OpenCOR include the incorporation of
GIT into OpenCOR to enable the upload of models to PMR, graphical
rendering of the model structure (using SVG), model building templates,
such as templates for creating Markov models, tools for parameter
estimation and tools for analysing model outputs.
The next release of CellML (1.2) will include the ability to specify a
probability distribution for a parameter value. Together with SED-ML,
this will allow OpenCOR to generate error bounds on the solutions,
corresponding to the specified parameter uncertainty.
Link to weblab.
These notes are currently being extended to include
- a discussion of system identification and parameter estimation
- more extensive discussion of membrane protein models
- CellML modules for signal transduction pathways
References
| [APJ15] | Garny A. and Hunter P.J. Opencor: a modular and interoperable approach to computational biology. Frontiers in Physiology, 2015. |
| [AYY] | Non A. Www.biomodels.org <http://www.biomodels.org>. YYYY. |
| [AAF52] | Hodgkin AL and Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. Journal of Physiology, 117:500–544, 1952. |
| [CPJ09] | Christie R. Nielsen P.M.F. Blackett S. Bradley C. and Hunter P.J. Fieldml: concepts and implementation. Philosophical Transactions of the Royal Society (London), A367(1895):1869–1884, 2009. |
| [CMEJ08] | Hunter PJ Cooling M and Crampin EJ. Modeling biological modularity with cellml. IET Systems Biology, 2:73–79, 2008. |
| [CJ15] | Waltemath D. Cooper J, Vik JO. A call for virtual experiments: accelerating the scientific process. Progress in Biophysics and Molecular Biolog, 117:99–106, 2015. |
| [D62] | Noble D. A modification of the hodgkin-huxley equations applicable to purkinje fibre action and pace-maker potentials. Journal of Physiology, 160:317–352, 1962. |
| [DPPJ03] | Cuellar A.A. Lloyd C.M. Nielsen P.F. Halstead M.D.B. Bullivant D.P. Nickerson D.P. and Hunter P.J. An overview of cellml 1.1, a biological model description language. SIMULATION: Transactions of the Society for Modeling and Simulation, 79(12):740–747, 2003. |
| [ea13] | Hunter P.J. et al. A vision and strategy for the virtual physiological human: 2012 update. Interface Focus, 2013. |
| [eal11] | Yu T. et al. The physiome model repository 2. Bioinformatics, 27:743–744, 2011. |
| [J97] | Wigglesworth J. Energy and life. Taylor & Francis Ltd, 1997. |
| [JH02] | Thompson JMT and Stewart HB. Nonlinear dynamics and chaos. Wiley, 2002. |
| [LCPF08] | Hunter PJ Lloyd CM, Lawson JR and Nielsen PF. The cellml model repository. Bioinformatics, 24:2122–2123, 2008. |
| [P13] | Britten R.D. Christie G.R. Little C. Miller A.K. Bradley C. Wu A. Yu T. Hunter P.J. Nielsen P. Fieldml, a proposed open standard for the physiome project for mathematical model representation. Med. Biol. Eng. Comput., 51(11):1191–1207, 2013. |
| [PJ04] | Hunter P.J. The iups physiome project: a framework for computational physiology. Progress in Biophysics and Molecular Biology, 85:551–569, 2004. |
| [VarYY] | Various. See www.cellml.org/about/publications for a more extensive list of publications on cellml and opencor. Various, YYYY. |
