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
Install and Launch OpenCOR
Download OpenCOR from www.opencor.ws.
Versions are available for Windows, OS X and Linux . Note that some aspects of this tutorial require
OpenCOR snapshot 2017-02-10 (or newer). Create a shortcut to the executable (found in the
bin directory) on your desktop and click on this to launch OpenCOR. A
window will appear that looks like Fig. 1(a).
OpenCOR application (a) Default positioning of dockable windows. (b) An
alternative configuration achieved by dragging and dropping the dockable
windows.
Dockable Windows
The central area is used to interact with files. By default, no files
are open, hence the OpenCOR logo is shown instead. To the sides, there
are dockable windows, which provide additional features. Those windows
can be dragged and dropped to the top or bottom of the central area as
shown in Figure 1(b) or they can be individually undocked or closed. All
closed panels can be re-displayed by enabling them in the View menu,
or by using the Tools menu Reset All option. The key combination
Control-spacebar removes (for less clutter) or restores these
two side panels .
Any of the subpanels (Physiome Model Repository, File Browser, and
File Organiser) can be closed with the top right delete button, and
then restored from the View .. Windows .. menu. Files can be dragged
and dropped into the File Organiser to create a local directory
structure for your files.
Plugins
OpenCOR has a plugin architecture and can be used with or
without a range of modules. These can be viewed under the Tools menu.
By default they are all included, as shown in Fig. 2. Information
about developing plugins for OpenCOR is also available.
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. 3.
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.
Simulation Panel
This area is used to set up the simulation settings.
Starting point - the value of the variable of integration (often
time) at which the simulation will begin. Leave this at 0.
Ending point - the point at which the simulation will end. Set to
100.
Point interval - the interval between data points on the variable of
integration. Set to 0.1.
Just above the Simulation panel are controls for running the
simulation. These are:
Run (
), Pause (
), Reset parameters (
),
Clear simulation data (
), Interval delay (
),
Add(
)/Subtract(
) graphical output windows
and Output solution to a CSV file (
).
For this model, we suggest that you create three graphical output
windows using the + button.
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. 4 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(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. 5 - 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. 5) 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. 7(a) for these parameters (a
decaying exponential) and in Fig. 7(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. 7 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. 9 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. 10. 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. 11).
Note
The simulation presented in Fig. 10 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
A model of ion channel gating and current: Introducing CellML units
A good example of a model based on a first order equation is the one
used by Hodgkin and Huxley [AAF52] to describe the gating behaviour of an
ion channel (see also next three sections). Before we describe the
gating behaviour of an ion channel, however, we need to explain the
concepts of the ‘Nernst potential’ and channel conductance.
An ion channel is a protein or protein complex embedded in the bilipid
membrane surrounding a cell and containing a pore through which an ion
\(Y^{+}\) (or \(Y^{-}\)) can pass when the channel is open. If
the concentration of this ion is
\(\left\lbrack Y^{+} \right\rbrack_{o}\) outside the cell and
\(\left\lbrack Y^{+} \right\rbrack_{i}\) inside the cell, the force
driving an ion through the pore is calculated from the change in
entropy.
Entropy \(S\) (\(J.K^{-1}\)) is a measure of the number of
microstates available to a system, as defined by Boltzmann’s equation
\(S = k_{B}\text{lnW}\), where \(W\) is the number of ways of
arranging a given distribution of microstates of a system and
\(k_{B}\) is Boltzmann’s constant. The driving force for ion
movement is the dispersal of energy into a more probable distribution
(see Fig. 12; cf. the second law of thermodynamics).
The energy change \(\Delta q\) associated with this change of
entropy \(\Delta S\) at temperature \(T\) is
\(\Delta q = T\Delta S\) (J).
For a given volume of fluid the number of microstates \(W\)
available to a solute (and hence the entropy of the solute) at a high
concentration is less than that for a low concentration. The
energy difference driving ion movement from a high ion concentration
\(\left\lbrack Y^{+} \right\rbrack_{i}\) (lower entropy) to a lower
ion concentration \(\left\lbrack Y^{+} \right\rbrack_{o}\) (higher
entropy) is therefore
\(\Delta q = T\Delta S = k_{B}T\left( \ln{\left\lbrack Y^{+} \right\rbrack_{o} - \ln\left\lbrack Y^{+} \right\rbrack_{i}} \right) = k_{B}T\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\)
(\(J.ion^{-1}\))
or
\(\Delta Q = RT\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\)
(\(J.mol^{-1}\)).
\(R = k_{B}N_{A} \approx 1.34x10^{-23} (J.K^{-1}) \text{x}
6.02x10^{23} (mol^{-1}) \approx 8.4 (J.mol^{-1}K^{-1})\)
is the ‘universal gas constant’.
At 25°C (298K), \(\text{RT} \approx 2.5 kJ.mol^{-1}\).
Every positively charged ion that crosses the membrane raises the
potential difference and produces an electrostatic driving force that
opposes the entropic force (see Fig. 13). To move an electron of
charge e (\(\approx 1.6x10^{-19}\) C) through a voltage change of
\(\Delta\phi\) (V) requires energy \(e\Delta\phi\) (J) and
therefore the energy needed to move an ion \(Y^{+}\) of valence
z=1 (the number of charges per ion) through a voltage change of
\(\Delta\phi\) is \(\text{ze}\Delta\phi\)
(\(J.ion^{-1}\)) or
\(\text{ze}N_{A}\Delta\phi\) (\(J.mol^{-1}\)). Using Faraday’s
constant \(F = eN_{A}\), where
\(F \approx 0.96x10^{5} C.mol^{-1}\), the change in energy
density at the macroscopic scale is \(\text{zF}\Delta\phi\)
(\(J.mol^{-1}\)).
No further movement of ions takes place when the force for entropy
driven ion movement exactly equals the opposing electrostatic driving
force associated with charge movement:
\(\text{zF}\Delta\phi = \text{RT}\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\)
(\(J.mol^{-1}\)) or
\(\Delta\phi = E_{Y} = \frac{\text{RT}}{\text{zF}}\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\)
(\(J.C^{-1}\) or V)
where \(E_{Y}\) is the ‘equilibrium’ or ‘Nernst’ potential for
\(Y^{+}\). At 25°C (298K),
\(\frac{\text{RT}}{F} = \frac{2.5x10^{3}\ }{0.96x10^{5}} (J.C^{-1}) \approx 25mV\).
For an open channel the electrochemical current flow is driven by the
open channel conductance \({\overset{\overline{}}{g}}_{Y}\) times
the difference between the transmembrane voltage \(V\) and the
Nernst potential for that ion:
\({\overset{\overline{}}{i}}_{Y}\mathbf{=}{\overset{\overline{}}{g}}_{Y}\left( V - E_{Y} \right)\).
This defines a linear current-voltage relation (‘Ohms law’) as shown in
Fig. 14. The gates to be discussed below modify this open channel
conductance.
To describe the time dependent transition between the closed and open
states of the channel, Hodgkin and Huxley introduced the idea of channel
gates that control the passage of ions through a membrane ion channel.
If the fraction of gates that are open is y, the fraction of gates
that are closed is \(1-y\), and a first order ODE can be used to describe
the transition between the two states (see Fig. 15):
\(\frac{\text{dy}}{\text{dt}} = \alpha_{y}\left( 1 - y \right) - \beta_{y}\text{.y}\)
where \(\alpha_{y}\)is the opening rate and \(\beta_{y}\) is
the closing rate.
The solution to this ODE is
\(y = \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}} + Ae^{- \left( \alpha_{y} + \beta_{y} \right)t}\)
The constant \(A\) can be interpreted as
\(A = y\left( 0 \right) - \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}}\)
as in the previous example and, with \(y\left( 0 \right) = 0\) (i.e.
all gates initially shut), the solution looks like Fig. 16(a).
The experimental data obtained by Hodgkin and Huxley for the squid axon,
however, indicated that the initial current flow began more slowly
(Fig. 16(b)) and they modelled this by assuming that the ion channel had
\(\gamma\) gates in series so that conduction would only occur when
all gates were at least partially open. Since \(y\) is the
probability of a gate being open, \(y^{\gamma}\) is the probability
of all \(\gamma\) gates being open (since they are assumed to be
independent) and the current through the channel is
\(i_{Y} = {\overset{\overline{}}{i}}_{Y}y^{\gamma} = y^{\gamma}{\overset{\overline{}}{g}}_{Y}\left( V - E_{Y} \right)\)
where
\({\overset{\overline{}}{i}}_{Y}{= \overset{\overline{}}{g}}_{Y}\left( V - E_{Y} \right)\)
is the steady state current through the open gate.
We can represent this in OpenCOR with a simple extension of the first
order ODE model, but in developing this model we will also demonstrate
the way in which CellML deals with units.
Note that the decision to deal with units in CellML, rather than just
ignoring them or insisting that all equations are represented in
dimensionless form, was made in order to be able to check the
physical consistency of all terms in each equation.
There are seven base physical quantities defined by the International d’Unités (SI). These are (with their SI units):
All other units are derived from these seven. Additional derived units
that CellML defines intrinsically (with their dependence on previously
defined units) are: Hz (\(s^{-1}\)); Newton, N
(\(kg.m.s^{-1}\)); Joule, J (\(N.m\)); Pascal, Pa (\(N.m^{-2}\));
Watt, W (\(J.s^{-1}\)); Volt, V (\(W.A^{-1}\)); Siemen, S
(\(A.V^{-1}\)); Ohm, \(\Omega\) (\(V.A^{-1}\)); Coulomb, C
(\(s.A\)); Farad, F (\(C.V^{-1}\)); Weber, Wb (\(V.s\)); and Henry,
H (\(Wb.A^{-1}\)). Multiples and fractions of these are defined as
follows:
|
Prefix |
|
deca |
hecto |
kilo |
mega |
giga |
tera |
peta |
exa |
zetta |
yotta |
Multiples |
Symbol |
|
da |
h |
k |
M |
G |
T |
P |
E |
Z |
Y |
|
Factor |
\(10^0\) |
\(10^{1}\) |
\(10^{2}\) |
\(10^{3}\) |
\(10^{6}\) |
\(10^{9}\) |
\(10^{12}\) |
\(10^{15}\) |
\(10^{18}\) |
\(10^{21}\) |
\(10^{24}\) |
|
Prefix |
|
deci |
centi |
milli |
micro |
nano |
pico |
femto |
atto |
zepto |
yocto |
Fractions |
Symbol |
|
d |
c |
m |
μ |
n |
p |
f |
a |
z |
y |
|
Factor |
\(10^{0}\) |
\(10^{-1}\) |
\(10^{-2}\) |
\(10^{-3}\) |
\(10^{-6}\) |
\(10^{-9}\) |
\(10^{-12}\) |
\(10^{-15}\) |
\(10^{-18}\) |
\(10^{-21}\) |
\(10^{-24}\) |
Units for this model, with multiples and fractions, are illustrated in
the following CellML Text code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32 | def model first_order_model as
def unit millisec as
unit second {pref: milli};
enddef;
def unit per_millisec as
unit second {pref: milli, expo: -1};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def unit milliS_per_cm2 as
unit siemens {pref: milli};
unit metre {pref: centi, expo: -2};
enddef;
def comp ion_channel as
var V: millivolt {init: 0};
var t: millisec {init: 0};
var y: dimensionless {init: 0};
var E_y: millivolt {init: -85};
var i_y: microA_per_cm2;
var g_y: milliS_per_cm2 {init: 36};
var gamma: dimensionless {init: 4};
var alpha_y: per_millisec {init: 1};
var beta_y: per_millisec {init: 2};
ode(y, t) = alpha_y*(1{dimensionless}-y)-beta_y*y;
i_y = g_y*pow(y, gamma)*(V-E_y);
enddef;
enddef;
|
Line 2: Define units for time as millisecs
Line 5: Define per_millisec units
Line 8: Define units for voltage as millivolts
Line 11: Define units for current as microAmps per cm2
Line 15: Define units for conductance as milliSiemens per cm2
Lines 20-28: Define units and initial conditions for variables
Line 29: Define ODE for gating variable y
Line 30: Define channel current
The solution of these equations for the parameters indicated above is
illustrated in Fig. 17.
The model of a gated ion channel presented here is used in the next two
sections for the neural potassium and sodium channels and then in
Section 11 for cardiac ion channels. The gates make the channel
conductance time dependent and, as we will see in the next section, the
experimentally observed voltage dependence of the gating rate constants
\(\alpha_{y}\) and \(\beta_{y}\) means that the channel
conductance (including the open channel conductance) is voltage
dependent. For a partially open channel (\(y < 1\)), the steady
state conductance is
\(\left( y_{\infty} \right)^{\gamma}{.\overset{\overline{}}{g}}_{Y}\),
where \(y_{\infty} = \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}}\).
Moreover the gating time constants
\(\tau = \frac{1}{\alpha_{y} + \beta_{y}}\) are therefore also
voltage dependent. Both of these voltage dependent factors of ion
channel gating are important in explaining channel properties, as we
show now for the neural potassium and sodium ion channels.
Footnotes
A model of the potassium channel: Introducing CellML components and connections
We now deal specifically with the application of the previous model to
the Hodgkin and Huxley (HH) potassium channel. Following the convention
introduced by Hodgkin and Huxley, the gating variable for the potassium
channel is \(n\) and the number of gates in series is
\(\gamma = 4\), therefore
\(i_{K} = \bar{i_K}n^{4} = n^{4}\bar{g}_{K}\left( V - E_{K} \right)\)
where \(\bar{g}_{K} = \ 36 \text{mS.cm}^{-2}\),
and with intra- and extra-cellular concentrations
\(\left\lbrack K^{+} \right\rbrack_{i} = 90\text{mM}\) and
\(\left\lbrack K^{+} \right\rbrack_{o} = 3\text{mM}\), respectively, the
Nernst potential for the potassium channel (\(z = 1\) since one +ve charge on
\(K^{+}\)) is
\(E_{k} = \frac{\text{RT}}{\text{zF}}\ ln\frac{\left\lbrack K^{+} \right\rbrack_{o}}{\left\lbrack K^{+} \right\rbrack_{i}} = 25\ ln\frac{3}{90} = - 85\text{mV}\).
As noted above, this is called the equilibrium potential since it is
the potential across the cell membrane when the channel is open but no
current is flowing because the electrostatic driving force from the
potential (voltage) difference between internal and external ion charges
is exactly matched by the entropic driving force from the ion
concentration difference. \(n^{4}\bar{g}_{K}\) is
the channel conductance.
The gating kinetics are described (as before) by
\(\frac{\text{dn}}{\text{dt}} = \alpha_{n}\left( 1 - n \right) - \beta_{n}\text{.n}\)
with time constant \(\tau_{n} = \frac{1}{\alpha_{n} + \beta_{n}}\)
(see A simple first order ODE).
The main difference from the gating model in our previous example is
that Hodgkin and Huxley found it necessary to make the rate constants
functions of the membrane potential \(V\) (see Fig. 18) as
follows:
\(\alpha_{n} = \frac{- 0.01\left( V + 65 \right)}{e^{\frac{- \left( V + 65 \right)}{10}} - 1}\);
\(\beta_{n} = 0.125e^{\frac{- \left( V + 75 \right)}{80}}\) .
Note that under steady state conditions when
\(t \rightarrow \infty\) and
\(\frac{\text{dn}}{\text{dt}} \rightarrow 0\),
\(\left. \ n \right|_{t = \infty} = n_{\infty} = \frac{\alpha_{n}}{\alpha_{n} + \beta_{n}}\).
The voltage dependence of the steady state channel conductance is then
\(g_{\text{SS}} = \left( \frac{\alpha_{n}}{\alpha_{n} + \beta_{n}} \right)^{4}.\bar{g}_{Y}\).
(see Fig. 18). The steady state current-voltage relation for the
channel is illustrated in Fig. 19.
These equations are captured with OpenCOR CellML Text view (together
with the previous unit definitions) below. But first we need
to explain some further CellML concepts.
We introduced CellML units above. We now need to introduce three
more CellML constructs: components, connections (mappings
between components) and groups. For completeness we also show one
other construct in Fig. 20, imports, that will be used later in A model of the nerve action potential: Introducing CellML imports.
Defining components serves two purposes: it preserves a modular
structure for CellML models, and allows these component modules to be
imported into other models, as we will illustrate later [DPPJ03]. For the
potassium channel model we define components representing (i) the
environment, (ii) the potassium channel conductivity, and (iii) the
dynamics of the n-gate.
Since certain variables (t, V and n) are shared between components, we
need to also define the component maps as indicated in the CellML Text
view below.
The CellML Text code for the potassium ion channel model is as
follows:
Potassium_ion_channel.cellml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80 | def model potassium_ion_channel as
def unit millisec as
unit second {pref: milli};
enddef;
def unit per_millisec as
unit second {pref: milli, expo: -1};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit per_millivolt as
unit millivolt {expo: -1};
enddef;
def unit per_millivolt_millisec as
unit per_millivolt;
unit per_millisec;
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def unit milliS_per_cm2 as
unit siemens {pref: milli};
unit metre {pref: centi, expo: -2};
enddef;
def unit mM as
unit mole {pref: milli};
enddef;
def comp environment as
var V: millivolt { pub: out};
var t: millisec {pub: out};
V = sel
case (t > 5 {millisec}) and (t < 15 {millisec}):
-85.0 {millivolt};
otherwise:
0.0 {millivolt};
endsel;
enddef;
def group as encapsulation for
comp potassium_channel incl
comp potassium_channel_n_gate;
endcomp;
enddef;
def comp 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 Ko: mM {init: 3};
var Ki: mM {init: 90};
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;
def comp potassium_channel_n_gate as
var V: millivolt {pub: in};
var t: millisec {pub: in};
var n: dimensionless {init: 0.325, pub: out};
var alpha_n: per_millisec;
var beta_n: per_millisec;
alpha_n = 0.01{per_millivolt_millisec}*(V+10{millivolt})
/(exp((V+10{millivolt})/10{millivolt})-1{dimensionless});
beta_n = 0.125{per_millisec}*exp(V/80{millivolt});
ode(n, t) = alpha_n*(1{dimensionless}-n)-beta_n*n;
enddef;
def map between environment and potassium_channel for
vars V and V;
vars t and t;
enddef;
def map between potassium_channel and
potassium_channel_n_gate for
vars V and V;
vars t and t;
vars n and n;
enddef;
enddef;
|
Lines 2-28: Define units.
Lines 29-38: Define component ‘environment’.
Lines 32-37: Define voltage step.
Lines 39-43: Define encapsulation of ‘n_gate’.
Lines 44-58: Define component ‘potassium_channel’.
Lines 59-69: Define component ‘potassium_channel_n_gate’.
Lines 70-79: Define mappings between components for variables that are shared between these components.
Note that several other features have been added:
the event control select case which indicates that the voltage is
specified to jump from 0 mV to -85 mV at t = 5 ms then back to 0 mV at
t = 15 ms. This is only used here in order to test the K channel model;
when the potassium_channel component is later imported into a neuron
model, the environment component is not imported.
the use of encapsulation to embed the
potassium_channel_n_gate inside the potassium_channel.
This avoids the need to establish mappings from environment to
potassium_channel_n_gate since the gate component is entirely
within the channel component.
the use of \(\left\{ pub:in \right\}\) and
\(\left\{ pub:out \right\}\) to indicate which variables are
either supplied as inputs to a component or produced as outputs from
a component. Any variables not labelled as in or out are
local variables or parameters defined and used only within that
component. Public (and private) interfaces are discussed in more
detail in the next section.
We now use OpenCOR, with Ending point 40 and Point interval 0.1, to
solve the equations for the potassium channel under a voltage step
condition in which the membrane voltage is clamped initially at 0mV and
then stepped down to -85mV for 10ms before being returned to 0mV. At
0mV, the steady state value of the n gate is
\(n_{\infty} = \frac{\alpha_{n}}{\alpha_{n} + \beta_{n}} =\) 0.324
and, at -85mV, \(n_{\infty} = \ \)0.945.
The voltage traces are shown at the top of Figure 21. The n-gate
response, shown next, is to open further from its partially open value
of \(n =\)0.324 at 0mV and then plateau at an almost fully open
state of \(n =\)0.945 at the Nernst potential -85mV before closing
again as the voltage is stepped back to 0mV. Note that the gate opening
behaviour (set by the voltage dependence of the \(\alpha_{n}\)
opening rate constant) is faster than the closing behaviour (set by the
voltage dependence of the \(\beta_{n}\) closing rate constant). The
channel conductance (\(= n^{4}\bar{g}_K\)) is
shown next – note the initial s-shaped conductance increase caused by
the \(n^{4}\) (four gates in series) effect on conductance. Finally
the channel current \(i_{K} =\) conductance x
\(\left( V - E_{K} \right)\) is shown at the bottom. Because the
voltage is clamped at the Nernst potential (-85mV) during the period
when the gate is opening, there is no current flow, but when the voltage
is stepped back to 0mV, the open gates begin to close and the
conductance declines but now there is a voltage gradient to drive an
outward (positive) current flow through the partially open channel –
albeit brief since the channel is closing.
Note that the CellML Text code above includes the Nernst equation with
its dependence on the concentrations
\(\left\lbrack K^{+} \right\rbrack_{i}\)= 90mM and
\(\left\lbrack K^{+} \right\rbrack_{o}\)= 3mM. Try raising the
external potassium concentration to
\(\left\lbrack K^{+} \right\rbrack_{o}\)= 10mM – you will then see
the Nernst potential increase from -85mV to -55mV and a negative
(inward) current flowing during the period when the membrane voltage is
clamped to -85mV. The cell is now in a ‘hyperpolarised’ state because
the potential is less than the equilibrium potential.
Note that you can change a model parameter such as
\(\left\lbrack K^{+} \right\rbrack_{o}\) either by changing the
value in the left hand Parameters window (which leaves the file
unchanged) or by editing the CellML Text code (which does change the
file when you save from CellML Text view – which you have to do to see
the effect of that change.
This potassium channel model will be used later, along with a sodium
channel model and a leakage channel model, to form the Hodgkin-Huxley
neuron model, where the membrane ion channel current flows are coupled
to the equations governing current flow along the axon to generate an
action potential.
Footnotes
A model of the sodium channel: Introducing CellML encapsulation and interfaces
The HH sodium channel has two types of gate, an \(m\) gate (of which
there are 3) that is initially closed (\(m = 0\)) before activating
and inactivating back to the closed state, and an \(h\) gate that is
initially open (\(h = 1\)) before activating and inactivating back
to the open state. The short period when both types of gate are open
allows a brief window current to pass through the channel. Therefore,
\[i_{\text{Na}} = \bar{i}_{\text{Na}}m^{3}h = m^{3}\text{h.}\bar{g}_{\text{Na}}\left( V - E_{\text{Na}} \right)\]
where \(\bar{g}_{\text{Na}} = \ \)120
mS.cm-2, and with
\(\left\lbrack \text{Na}^{+} \right\rbrack_{i}\)= 30mM and
\(\left\lbrack \text{Na}^{+} \right\rbrack_{o}\)= 140mM, the
Nernst potential for the sodium channel (z=1) is
\[E_{\text{Na}} = \frac{\text{RT}}{\text{zF}}ln\frac{\left\lbrack \text{Na}^{+} \right\rbrack_{o}}{\left\lbrack \text{Na}^{+} \right\rbrack_{i}} = 25\ ln\frac{140}{30} = 35\text{mV}.\]
The gating kinetics are described by
\[\frac{\text{dm}}{\text{dt}} = \alpha_{m}\left( 1 - m \right) - \beta_{m}\text{.m};
\frac{\text{dh}}{\text{dt}} = \alpha_{h}\left( 1 - h \right) - \beta_{h}\text{.h}\]
where the voltage dependence of these four rate constants is determined
experimentally to be
\[\alpha_{m} = \frac{- 0.1\left( V + 50 \right)}{e^{\frac{- \left( V + 50 \right)}{10}} - 1};
\beta_{m} = 4e^{\frac{- \left( V + 75 \right)}{18}};
\alpha_{h} = 0.07e^{\frac{- \left( V + 75 \right)}{20}};
\beta_{h} = \frac{1}{e^{\frac{- \left( V + 45 \right)}{10}} + 1}.\]
Before we construct a CellML model of the sodium channel, we first
introduce some further CellML concepts that help deal with the
complexity of biological models: first the use of encapsulation groups
and public and private interfaces to control the visibility of
information in modular CellML components. To understand encapsulation,
it is useful to use the terms ‘parent’, ‘child’ and ‘sibling’.
def group as encapsulation for
comp sodium_channel incl
comp sodium_channel_m_gate;
comp sodium_channel_h_gate;
endcomp;
enddef;
We define the CellML components sodium_channel_m_gate and
sodium_channel_h_gate below. Each of these components has its own
equations (voltage-dependent gates and first order gate kinetics) but
they are both parts of one protein – the sodium channel – and it is
useful to group them into one sodium_channel component as shown above:
We can then talk about the sodium channel as the parent of two children:
the m gate and the h gate, which are therefore siblings. A private
interface allows a parent to talk to its children and a public
interface allows siblings to talk among themselves and to their parents
(see Fig. 22).
The OpenCOR CellML Text for the HH sodium ion channel is given below.
Sodium_ion_channel.cellml
def model sodium_ion_channel as
def unit millisec as
unit second {pref: milli};
enddef;
def unit per_millisec as
unit second {pref: milli, expo: -1};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit per_millivolt as
unit millivolt {expo: -1};
enddef;
def unit per_millivolt_millisec as
unit per_millivolt;
unit per_millisec;
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def unit milliS_per_cm2 as
unit siemens {pref: milli};
unit metre {pref: centi, expo: -2};
enddef;
def unit mM as
unit mole {pref: milli};
enddef;
def comp environment as
var V: millivolt {pub: out};
var t: millisec {pub: out};
V = sel
case (t > 5 {millisec}) and (t < 15 {millisec}):
-20.0 {millivolt};
otherwise:
-85.0 {millivolt};
endsel;
enddef;
def group as encapsulation for
comp sodium_channel incl
comp sodium_channel_m_gate;
comp sodium_channel_h_gate;
endcomp;
enddef;
def comp sodium_channel as
var V: millivolt {pub: in, priv: out};
var t: millisec {pub: in, priv: out };
var m: dimensionless {priv: in};
var h: dimensionless {priv: in};
var g_Na: milliS_per_cm2 {init: 120};
var i_Na: microA_per_cm2 {pub: out};
var Nao: mM {init: 140};
var Nai: mM {init: 30};
var RTF: millivolt {init: 25};
var E_Na: millivolt;
var Na_conductance: milliS_per_cm2 {pub: out};
E_Na=RTF*ln(Nao/Nai);
Na_conductance = g_Na*pow(m, 3{dimensionless})*h;
i_Na= Na_conductance*(V-E_Na);
enddef;
def comp sodium_channel_m_gate as
var V: millivolt {pub: in};
var t: millisec {pub: in};
var alpha_m: per_millisec;
var beta_m: per_millisec;
var m: dimensionless {init: 0.05, pub: out};
alpha_m = 0.1{per_millivolt_millisec}*(V+25{millivolt})
/(exp((V+25{millivolt})/10{millivolt})-1{dimensionless});
beta_m = 4{per_millisec}*exp(V/18{millivolt});
ode(m, t) = alpha_m*(1{dimensionless}-m)-beta_m*m;
enddef;
def comp sodium_channel_h_gate as
var V: millivolt {pub: in};
var t: millisec {pub: in};
var alpha_h: per_millisec;
var beta_h: per_millisec;
var h: dimensionless {init: 0.6, pub: out};
alpha_h = 0.07{per_millisec}*exp(V/20{millivolt});
beta_h = 1{per_millisec}/(exp((V+30{millivolt})/10{millivolt})+1{dimensionless});
ode(h, t) = alpha_h*(1{dimensionless}-h)-beta_h*h;
enddef;
def map between environment and sodium_channel for
vars V and V;
vars t and t;
enddef;
def map between sodium_channel and sodium_channel_m_gate for
vars V and V;
vars t and t;
vars m and m;
enddef;
def map between sodium_channel and sodium_channel_h_gate for
vars V and V;
vars t and t;
vars h and h;
enddef;
enddef;
The results of the OpenCOR computation, with Ending point 40 and
Point interval 0.1, are shown in Fig. 23 with plots \(V\left( t \right)\), \(m\left( t \right)\),
\(h\left( t \right)\), \(g_{\text{Na}}\left( t \right)\) and
\(i_{\text{Na}}(t)\) for voltage steps from (a) -85mV to -20mV, (b) -85mV to 0mV and (c) -85mV to 20mV. There are several
things to note:
The kinetics of the m-gate are much faster than the h-gate.
The opening behaviour is faster as the voltage is stepped to higher
values since \(\tau = \frac{1}{\alpha_{n} + \beta_{n}}\)
reduces with increasing V (see Fig. 18).
The sodium channel conductance rises (activates) and then falls
(inactivates) under a positive voltage step from rest since the
three m-gates turn on but the h-gate turns off and the conductance
is a product of these. Compare this with the potassium channel
conductance shown in Fig. 21 which is only reduced back to zero
by stepping the voltage back to its resting value – i.e.
deactivating it.
The only time current \(i_{\text{Na}}\) flows through the
sodium channel is during the brief period when the m-gate is
rapidly opening and the much slower h-gate is beginning to close. A
small current flows during the reverse voltage step but this is at
a time when the h-gate is now firmly off so the magnitude is very
small.
The large sodium current \(i_{\text{Na}}\) is an inward current
and hence negative.
Note that the bottom trace does not quite line up at t=0 because the
values shown on the axes are computed automatically and hence can take
more or less space depending on their magnitude.
Footnotes
A model of the nerve action potential: Introducing CellML imports
Here we describe the first (and most famous) model of nerve fibre
electrophysiology based on the membrane ion channels that we have
discussed in the last two sections. This is the work by Alan Hodgkin and
Andrew Huxley in 1952 [AAF52] that won them (together with John Eccles) the
1963 Noble prize in Physiology or Medicine for “their discoveries
concerning the ionic mechanisms involved in excitation and inhibition in
the peripheral and central portions of the nerve cell membrane”.
Cable equation
The cable equation was developed in 1890 to predict the
degradation of an electrical signal passing along the transatlantic
cable. It is derived as follows:
If the voltage is raised at the left hand end of the cable (shown by the
deep red in Fig. 24), a current \(i_{a}\) (A) will flow that
depends on the voltage gradient, given by
\(\frac{\partial V}{\partial x}\) (\(V.m^{-1}\)) and the resistance
\(r_{a}\) (\(\Omega.m^{-1}\)), Ohm’s law gives
\(- \frac{\partial V}{\partial x} = r_{a}i_{a}\) . But if the cable
leaks current \(i_{m}\) (\(A.m^{-1}\)) per unit length of cable,
conservation of current gives
\(\frac{\partial i_{a}}{\partial x} = i_{m}\) and therefore,
substituting for \(i_{a}\) ,
\(\frac{\partial}{\partial x}\left( - \frac{1}{r_{a}}\frac{\partial V}{\partial x} \right) = i_{m}\)
. There are two sources of membrane current \(i_{m}\) , one
associated with the capacitance \(C_{m}\)
(\(\approx 1\mu F/\text{cm}^{2}\)) of the membrane,
\(C_{m}\frac{\partial V}{\partial t}\), and one associated with
holes or channels in the membrane, \(i_{\text{leak}}\). Inserting
these into the RHS gives
\[\frac{\partial}{\partial x}\left( - \frac{1}{r_{a}}\frac{\partial V}{\partial x} \right) = i_{m} = C_{m}\frac{\partial V}{\partial t} + i_{\text{leak}}\]
Rearranging gives the cable equation (for constant \(r_{a}\)):
\[C_{m}\frac{\partial V}{\partial t} = - \frac{1}{r_{a}}\frac{\partial^{2}V}{\partial x^{2}} - i_{\text{leak}}\]
where all terms represent current density (current per membrane area)
and have units of \(\mu A/\text{cm}^{2}\).
Action potentials
The cable equation can be used to model the propagation of an action
potential along a neuron or any other excitable cell. The ‘leak’ current
is associated primarily with the inward movement of sodium ions through
the membrane ‘sodium channel’, giving the inward membrane current
\(i_{\text{Na}}\), and the outward movement of potassium ions
through a membrane ‘potassium channel’, giving the outward current
\(i_{K}\) (see Fig. 25). A further small leak current
\(i_{L} = g_{L}\left( V - E_{L} \right)\) associated with chloride
and other ions is also included.
When the membrane potential \(V\) rises due to axial current flow,
the Na channels open and the K channels close, such that the membrane
potential moves towards the Nernst potential for sodium. The subsequent
decline of the Na channel conductance and the increasing K channel
conductance as the voltage drops rapidly repolarises the membrane to its
resting potential of -85mV (see Fig. 26).
We can neglect the term
(\(- \frac{1}{r_{a}}\frac{\partial^{2}V}{\partial x^{2}}\)) (the
rate of change of axial current along the cable) for the present models
since we assume the whole cell is clamped with an axially uniform
potential. We can therefore obtain the membrane potential \(V\) by
integrating the first order ODE
\[\frac{\text{dV}}{\text{dt}} = - \left( i_{\text{Na}} + \ i_{K} + i_{L} \right)/C_{m}.\]
We use this example to demonstrate the importing feature of CellML.
CellML imports are used to bring a previously defined CellML model of
a component into the new model (in this case the Na and K channel
components defined in the previous two sections, together with a leakage
ion channel model specified below). Note that importing a component
brings the children components with it along with their connections and
units, but it does not bring the siblings of that component with it.
To establish a CellML model of the HH equations we first lay out the
model components with their public and private interfaces (Fig. 28).
The HH model is the top level model. The CellML Text code for the HH
model, together with the leakage_channel model, is given below. The imported potassium_ion_channel model and
sodium_ion_channel model are unchanged from the previous sections
HH.cellml
def model HH as
def import using "sodium_ion_channel.cellml" for
comp Na_channel using comp sodium_channel;
enddef;
def import using "potassium_ion_channel.cellml" for
comp K_channel using comp potassium_channel;
enddef;
def import using "leakage_ion_channel.cellml" for
comp L_channel using comp leakage_channel;
enddef;
def unit millisec as
unit second {pref: milli};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def unit microF_per_cm2 as
unit farad {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def group as encapsulation for
comp membrane incl
comp Na_channel;
comp K_channel;
comp L_channel;
endcomp;
enddef;
def comp environment as
var V: millivolt {init: -85, pub: out};
var t: millisec {pub: out};
enddef;
def map between environment and membrane for
vars V and V;
vars t and t;
enddef;
def map between membrane and Na_channel for
vars V and V;
vars t and t;
vars i_Na and i_Na;
enddef;
def map between membrane and K_channel for
vars V and V;
vars t and t;
vars i_K and i_K;
enddef;
def map between membrane and L_channel for
vars V and V;
vars i_L and i_L;
enddef;
def comp membrane as
var V: millivolt {pub: in, priv: out};
var t: millisec {pub: in, priv: out};
var i_Na: microA_per_cm2 {pub: out, priv: in};
var i_K: microA_per_cm2 {pub: out, priv: in};
var i_L: microA_per_cm2 {pub: out, priv: in};
var Cm: microF_per_cm2 {init: 1};
var i_Stim: microA_per_cm2;
var i_Tot: microA_per_cm2;
i_Stim = sel
case (t >= 1{millisec}) and (t <= 1.2{millisec}):
100{microA_per_cm2};
otherwise:
0{microA_per_cm2};
endsel;
i_Tot = i_Stim + i_Na + i_K + i_L;
ode(V,t) = -i_Tot/Cm;
enddef;
enddef;
Leakage_ion_channel
def model leakage_ion_channel as
def unit millisec as
unit second {pref: milli};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit per_millivolt as
unit millivolt {expo: -1};
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: -2};
enddef;
def unit milliS_per_cm2 as
unit siemens {pref: milli};
unit metre {pref: centi, expo: -2};
enddef;
def comp environment as
var V: millivolt {init: 0, pub: out};
var t: millisec {pub: out};
enddef;
def map between leakage_channel and environment for
vars V and V;
enddef;
def comp leakage_channel as
var V: millivolt {pub: in};
var i_L: microA_per_cm2 {pub: out};
var g_L: milliS_per_cm2 {init: 0.3};
var E_L: millivolt {init: -54.4};
i_L = g_L*(V-E_L);
enddef;
enddef;
Note that the CellML Text code for the potassium channel is Potassium_ion_channel.cellml
and for the sodium channel is Sodium_ion_channel.cellml.
Note that the only units that need to be defined for this top level HH
model are the ones explicitly required for the membrane component. All
the other units, required for the various imported sub-models, are
imported along with the imported components.
The results generated by the HH model are shown in Fig. 29.
Important note
It is often convenient to have the sub-models – in this case the
sodium_ion_channel.cellml model, the potassium_ion_channel.cellml
model and the leakage_ion_channel.cellml model - loaded into OpenCOR
at the same time as the high level model (HH.cellml), as shown in Fig. 30
. If you make changes to a model in the CellML Text view, you must
save the file (CTRL-S) before running a new simulation since the
simulator works with the saved model. Furthermore, a change to a
sub-model will only affect the high level model which imports it if you
also save the high level model (or use the Reload option under the
File menu). An asterisk appears next to the name of a file when a change
has been made and the file has not been saved. The asterisk disappears
when the file is saved.
Footnotes
A model of the cardiac action potential: Importing units and parameters
We now examine the Noble 1962 model [D62] that applied the Hodgkin-Huxley
approach to cardiac cells and thereby initiated the development of a
long line of cardiac cell models that, in their human cell formulation,
are now used clinically and are the most sophisticated models of any
cell type. It was the incorporation of these models into whole heart
bioengineering models that initiated the Physiome Project. We also
illustrate the use of imported units and imported parameter sets.
Cardiac cells have similar gradients of potassium and sodium ions that
operate in a similar way to neurons (as do all electrically active
cells). There is one major difference, however, in the potassium channel
that holds the cells in their resting state at -85mV (HH neuron) or
-100mV (cardiac Purkinje cells). This difference is illustrated in
Fig. 31(a). When the membrane potential is raised above the equilibrium
potential for potassium, the cardiac channel conductance shown by the
dashed line drops to nearly zero – i.e. it is an inward rectifier
since it rectifies (‘cuts off’) the outward current that otherwise would
have flowed through the channel at that potential. This is an
evolutionary adaptation of the potassium channel to avoid loss of
potassium ions out of the cell during the long plateau phase of the
cardiac action potential (Fig. 31(b)) needed to give the heart time to
contract. This evolutionary change saves the additional energy that
would otherwise be needed to pump potassium ions back into the cell, but
this Faustian “pact with the devil” is also the reason the heart is so
susceptible to conduction failure (more on this later). To explain his
data on Purkinje cells Noble [D62] postulated the existence of two inward
rectifier potassium channels, one with a conductance \(g_{K1}\) that
showed voltage dependence but no significant time dependence and another
with conductance \(g_{K2}\) that showed less severe rectification
with time dependent gating similar to the HH four-gated potassium
channel.
To model the cardiac action potential in Purkinje fibres (a cardiac cell
specialised for rapid conduction from the atrio-ventricular node to the
apical ventricular myocardial tissue), Noble [D62] proposed two potassium
channels (one of these being the inwardly rectifying potassium channel
described above and the other called the delayed potassium channel), one
sodium channel (very similar to the HH neuronal sodium channel) and one
leakage channel (also similar to the HH one).
The equations for these are as follows: (as for the HH model, time is in
ms, voltages are in mV, concentrations are in mM, conductances are in
mS, currents are in µA and capacitance is in µF).
Inward rectifying \(\mathbf{i}_{\mathbf{K}\mathbf{1}}\)
potassium channel (voltage dependent only)
\[\begin{split}i_{K1} &=\ g_{K1}\left( V - E_{K} \right),\text{ with }
E_{K} = \frac{\text{RT}}{\text{zF}}ln\frac{\left\lbrack K^{+} \right\rbrack_{o}}{\left\lbrack K^{+} \right\rbrack_{i}} = 25ln\frac{2.5}{140} = - 100\text{mV}. \\
g_{K1} &=\ 1.2e^{\frac{- \left( V + 90 \right)}{50}} + 0.015e^{\frac{\left( V + 90 \right)}{60}}\end{split}\]
Inward rectifying \(\mathbf{i}_{\mathbf{K}\mathbf{2}}\)
potassium channel (voltage and time dependent)
\[\begin{split}i_{K2} &=\ g_{K2}\left( V - E_{K} \right) \\
g_{K2} &=\ 1.2n^{4} \\
\frac{\text{dn}}{\text{dt}} &=\ \alpha_{n}\left( 1 - n \right) - \beta_{n}\text{.n},
\text{ where } \alpha_{n} = \frac{- 0.0001\left( V + 50 \right)}{e^{\frac{- \left( V + 50 \right)}{10}} - 1} \text{ and } \beta_{n} = 0.002e^{\frac{- \left( V + 90 \right)}{80}}.\end{split}\]
Note that the rate constants here reflect a much slower onset of the
time dependent change in conductance than in the HH potassium channel.
Sodium channel
\[\begin{split}i_{\text{Na}} &=\ \left( g_{\text{Na}} + 140 \right)\left( V - E_{\text{Na}} \right), \text{ with } E_{\text{Na}} = \frac{\text{RT}}{\text{zF}}ln\frac{\left\lbrack \text{Na}^{+} \right\rbrack_{o}}{\left\lbrack \text{Na}^{+} \right\rbrack_{i}} = 25ln\frac{140}{30} = 35\text{mV}. \\
g_{\text{Na}} &=\ m^{3}\text{h.}g_{Na\_ max} \text{ where } g_{Na\_ max} = 400\text{mS}. \\
\frac{\text{dm}}{\text{dt}} &=\ \alpha_{m}\left( 1 - m \right) - \beta_{m}\text{.m}, \text{ where } \alpha_{m} = \frac{- 0.1\left( V + 48 \right)}{e^{\frac{- \left( V + 48 \right)}{15}} - 1} \text{ and } \beta_{m} = \frac{0.12\left( V + 8 \right)}{e^{\frac{\left( V + 8 \right)}{5}} - 1} \\
\frac{\text{dh}}{\text{dt}} &=\ \alpha_{h}\left( 1 - h \right) - \beta_{h}\text{.h},
\text{ where } \alpha_{h} = 0.17e^{\frac{- \left( V + 90 \right)}{20}}\text{ and } \beta_{h} = \frac{1}{1 + e^{\frac{- \left( V + 42 \right)}{10}}}\end{split}\]
Leakage channel
\[i_{\text{leak}} = g_{L}\left( V - E_{L} \right), \text{ with }
E_{L} = - 60mV \text{ and } g_{L} = 0.075\text{mS}.\]
Membrane equation
\[\frac{\text{dV}}{\text{dt}} = - \left( i_{\text{Na}} + i_{K1} + i_{K2} + i_{\text{leak}} \right)/C_{m}\text{ where } C_{m} = 12\mu\text{F}.\]
Fig. 32 shows the structure of the model, including separate files for
units, parameters, and the three ion channels (the two potassium
channels are lumped together). We include the Nernst equations
dependence on potassium and sodium ion concentrations in order to
demonstrate the use of parameter values, defined in a separate
parameters file, that are read in at the top (whole cell model) level
and passed down to the individual ion channel models.
The CellML Text code for all six files is shown on the following two
pages. The arrows indicate the imports (appropriately colour coded for
units, components, and parameters).
Graphical outputs from solution of the Noble 1962 model with OpenCOR for
5000ms are shown in Fig. 32. Interpretation of the model outputs is given in the Fig. 32 legend.
The Noble62 model was developed further by Noble and others to include
additional sodium and potassium channels, calcium channels (needed for
excitation-contraction coupling), chloride channels and various ion
exchange mechanisms (Na/Ca, Na/H), co-transporters (Na/Cl, K/Cl) and
energy (ATP)-dependent pumps (Na/K, Ca) needed to model the observed
beat by beat changes in intracellular ion concentrations. These are
discussed further in Section 15.
Note
The downloadable links below are links to the raw text file that may be used for copying and pasting into OpenCOR. For the underlying CellML file that is suitable for opening with OpenCOR from disk obtain the xml file.
Raw text: Noble_1962.txt, XML file: Noble_1962.cellml.
def model Noble_1962 as
def import using "Noble62_Na_channel.xml" for
comp Na_channel using comp sodium_channel;
enddef;
def import using "Noble62_K_channel.xml" for
comp K_channel using comp potassium_channel;
enddef;
def import using "Noble62_L_channel.xml" for
comp L_channel using comp leakage_channel;
enddef;
def import using "Noble62_units.xml" for
unit mV using unit mV;
unit ms using unit ms;
unit nanoF using unit nanoF;
unit nanoA using unit nanoA;
enddef;
def import using "Noble62_parameters.xml" for
comp parameters using comp parameters;
enddef;
def map between parameters and membrane for
vars Ki and Ki;
vars Ko and Ko;
vars Nai and Nai;
vars Nao and Nao;
enddef;
def comp environment as
var t: ms {init: 0, pub: out};
enddef;
def group as encapsulation for
comp membrane incl
comp Na_channel;
comp K_channel;
comp L_channel;
endcomp;
enddef;
def comp membrane as
var V: mV {init: -85, pub: out, priv: out};
var t: ms {pub: in, priv: out};
var Cm: nanoF {init: 12000};
var Ki: mM {pub: in, priv: out};
var Ko: mM {pub: in, priv: out};
var Nai: mM {pub: in, priv: out};
var Nao: mM {pub: in, priv: out};
var i_Na: nanoA {pub: out, priv: in};
var i_K: nanoA {pub: out, priv: in};
var i_L: nanoA {pub: out, priv: in};
ode(V, t) = -(i_Na+i_K+i_L)/Cm;
enddef;
def map between environment and membrane for
vars t and t;
enddef;
def map between membrane and Na_channel for
vars V and V;
vars t and t;
vars Nai and Nai;
vars Nao and Nao;
vars i_Na and i_Na;
enddef;
def map between membrane and K_channel for
vars V and V;
vars t and t;
vars Ki and Ki;
vars Ko and Ko;
vars i_K and i_K;
enddef;
def map between membrane and L_channel for
vars V and V;
vars i_L and i_L;
enddef;
enddef;
Raw text: Noble62_units.txt, XML file Noble62_units.cellml.
def model Noble62_units as
def unit ms as
unit second {pref: milli};
enddef;
def unit per_ms as
unit second {pref: milli, expo: -1};
enddef;
def unit mV as
unit volt {pref: milli};
enddef;
def unit mM as
unit mole {pref: milli};
enddef;
def unit per_mV as
unit volt {pref: milli, expo: -1};
enddef;
def unit per_mV_ms as
unit mV {expo: -1};
unit ms {expo: -1};
enddef;
def unit microS as
unit siemens {pref: micro};
enddef;
def unit nanoF as
unit farad {pref: nano};
enddef;
def unit nanoA as
unit ampere {pref: nano};
enddef;
enddef;
Raw text: Noble62_parameters.txt, XML file Noble62_parameters.cellml.
def model Noble62_parameters as
def import using "Noble62_units.xml" for
unit mM using unit mM;
enddef;
def comp parameters as
var Ki: mM {init: 140, pub: out};
var Ko: mM {init: 2.5, pub: out};
var Nai: mM {init: 30, pub: out};
var Nao: mM {init: 140, pub: out};
enddef;
enddef;
Raw text: Noble62_Na_channel.txt, XML file Noble62_Na_channel.cellml.
def model sodium_ion_channel as
def import using "Noble62_units.xml" for
unit mV using unit mV;
unit ms using unit ms;
unit mM using unit mM;
unit per_ms using unit per_ms;
unit per_mV using unit per_mV;
unit per_mV_ms using unit per_mV_ms;
unit microS using unit microS;
unit nanoA using unit nanoA;
enddef;
def group as encapsulation for
comp sodium_channel incl
comp sodium_channel_m_gate;
comp sodium_channel_h_gate;
endcomp;
enddef;
def comp sodium_channel as
var V: mV {pub: in, priv: out};
var t: ms {pub: in, priv: out};
var g_Na_max: microS {init: 400000};
var g_Na: microS;
var E_Na: mV;
var m: dimensionless {priv: in};
var h: dimensionless {priv: in};
var Nai: mM {pub: in};
var Nao: mM {pub: in};
var RTF: mV {init: 25};
var i_Na: nanoA {pub: out};
E_Na = RTF*ln(Nao/Nai);
g_Na = pow(m, 3{dimensionless})*h*g_Na_max;
i_Na = (g_Na+140{microS})*(V-E_Na);
enddef;
def comp sodium_channel_m_gate as
var V: mV {pub: in};
var t: ms {pub: in};
var m: dimensionless {init: 0.01, pub: out};
var alpha_m: per_ms;
var beta_m: per_ms;
alpha_m = -0.10{per_mV_ms}*(V+48{mV})
/(exp(-(V+48{mV})/15{mV})-1{dimensionless});
beta_m = 0.12{per_mV_ms}*(V+8{mV})
/(exp((V+8{mV})/5{mV})-1{dimensionless});
ode(m, t)=alpha_m*(1{dimensionless}-m)-beta_m*m;
enddef;
def comp sodium_channel_h_gate as
var V: mV {pub: in};
var t: ms {pub: in};
var h: dimensionless {init: 0.8, pub: out};
var alpha_h: per_ms;
var beta_h: per_ms;
alpha_h = 0.17{per_ms}*exp(-(V+90{mV})/20{mV});
beta_h = 1.00{per_ms}
/(1{dimensionless}+exp(-(V+42{mV})/10{mV}));
ode(h, t) = alpha_h*(1{dimensionless}-h)-beta_h*h;
enddef;
def map between sodium_channel
and sodium_channel_m_gate for
vars V and V;
vars t and t;
vars m and m;
enddef;
def map between sodium_channel
and sodium_channel_h_gate for
vars V and V;
vars t and t;
vars h and h;
enddef;
enddef;
Raw text: Noble62_K_channel.txt, XML file Noble62_K_channel.cellml.
def model potassium_ion_channel as
def import using "Noble62_units.xml" for
unit mV using unit mV;
unit ms using unit ms;
unit mM using unit mM;
unit per_ms using unit per_ms;
unit per_mV using unit per_mV;
unit per_mV_ms using unit per_mV_ms;
unit microS using unit microS;
unit nanoA using unit nanoA;
enddef;
def group as encapsulation for
comp potassium_channel incl
comp potassium_channel_n_gate;
endcomp;
enddef;
def comp potassium_channel as
var V: mV {pub: in, priv: out};
var t: ms {pub: in, priv: out};
var n: dimensionless {priv: in};
var Ki: mM {pub: in};
var Ko: mM {pub: in};
var RTF: mV {init: 25};
var E_K: mV;
var g_K1: microS;
var g_K2: microS;
var i_K: nanoA {pub: out};
E_K = RTF*ln(Ko/Ki);
g_K1 = 1200{microS}*exp(-(V+90{mV})/50{mV})
+15{microS}*exp((V+90{mV})/60{mV});
g_K2 = 1200{microS}*pow(n, 4{dimensionless});
i_K = (g_K1+g_K2)*(V-E_K);
enddef;
def comp potassium_channel_n_gate as
var V: mV {pub: in};
var t: ms {pub: in};
var n: dimensionless {init: 0.01, pub: out};
var alpha_n: per_ms;
var beta_n: per_ms;
alpha_n = -0.0001{per_mV_ms}*(V+50{mV})
/(exp(-(V+50{mV})/10{mV})-1{dimensionless});
beta_n = 0.0020{per_ms}*exp(-(V+90{mV})/80{mV});
ode(n,t)= alpha_n*(1{dimensionless}-n)-beta_n*n;
enddef;
def map between environment
and potassium_channel for
vars V and V;
vars t and t;
enddef;
def map between potassium_channel and
potassium_channel_n_gate for
vars V and V;
vars t and t;
vars n and n;
enddef;
enddef;
Raw text: Noble62_L_channel.txt, XML file Noble62_L_channel.cellml.
def model leakage_ion_channel as
def import using "Noble62_units.xml" for
unit mV using unit mV;
unit ms using unit ms;
unit microS using unit microS;
unit nanoA using unit nanoA;
enddef;
def comp leakage_channel as
var V: mV {pub: in};
var g_L: microS {init: 75};
var E_L: mV {init: -60};
var i_L: nanoA {pub: out};
i_L = g_L*(V-E_L);
enddef;
enddef;
We have now covered all existing features of CellML and OpenCOR. But,
most importantly, you have learned ‘best practice’ for building CellML
models, including encapsulation of sub-components and a modular approach
in which units, parameters and model components are defined in separate
files that are imported into a composite model.
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. 34). 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. 35 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. 36. 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. 37).
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. 36. 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. 38, 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
The Physiome Model Repository and the link to bioinformatics
The Physiome Model Repository (PMR) [LCPF08] is the main online repository
for the IUPS Physiome Project, providing version and access controlled
repositories, called workspaces, for users to store their data.
Currently there are over 700 public workspaces and many
private workspaces in the repository. PMR also provides a mechanism to
create persistent access to specific revisions of a workspace, termed
exposures. Exposure plugins are available for specific types of data
(e.g. CellML or FieldML documents) which enable customizable views of
the data when browsing the repository via a web browser, or an
application accessing the repository’s content via web services.
More complete documentation describing how to use PMR is available in the PMR documentation: https://models.physiomeproject.org/docs.
The CellML models on models.physiomeproject.org are listed under 20 categories, shown below:
(numbers of exposures in each category are given besides the bar graph, correct as at early 2016)
Browse by category
Calcium Dynamics |
 140
|
Cardiovascular Circulation |
60 |
Cell Cycle |
38 |
Cell Migration |
2 |
Circadian Rhythms |
22 |
Electrophysiology |
230 |
Endocrine |
60 |
Excitation-Contraction Coupling |
22 |
Gene Regulation |
12 |
Hepatology |
29 |
Immunology |
55 |
Ion transport |
13 |
Mechanical Constitutive Laws |
19 |
Metabolism |
86 |
Myofilament Mechanics |
22 |
Neurobiology |
33 |
pH regulation |
2 |
PKPD |
11 |
Signal Transduction |
120 |
Synthetic Biology |
6 |
Note that searching of models can be done anywhere on the site using the
search box on the upper right hand corner. An important benefit of
ensuring that the models on the PMR are annotated is that models can
then be retrieved by a web-search using any of the annotated terms in
the models.
To illustrate the features of PMR, click on the Hund, Rudy 2004 (Basic)
model in the alphabetic listing of models under the Electrophysiology category.
The section labelled ‘Model Structure’ contains the journal paper
abstract and often a diagram of the model. This is shown for the
Hund-Rudy 2004 model in Fig. 40. This model, with over 22 separate
protein model components, is also a good example of why it is important
to build models from modular components [CMEJ08], and in particular the
individual ion channels for electrophysiology models.
There is a list of ‘Views Available’ for the CellML model on the
right hand side of the exposure page. The function of each of these
views is as follows:
Views Available
Documentation - Takes you to the main exposure page.
Model Metadata - Lists metadata including authors, title, journal,
Pubmed ID and model annotations.
Model Curation - Provides the curation status of the model. Note:
this is soon to be updated.
Mathematics - Displays all the mathematical equations contained in
the model.
Generated Code - Various codes (C, C-IDA, F77, MATLAB or Python) generated from
the model.
Cite this model - Provides details on how to cite use of the CellML
model.
Source view - Gives a full listing of the XML code for the model.
Launch with OpenCOR - Opens the model (or simulation experiment) in OpenCOR.
Note that CellML models are available under a Creative Commons
Attribution 3.0 Unported License. This means that you are free to:
Share — copy and redistribute the material in any medium or format
Adapt — remix, transform, and build upon the material
for any purpose, including commercial use.
The next stage of content development for PMR is to provide a list of
the modular components of these models each with their own exposure. For
example, models for each of the individual ion channels used in the
publication-based electrophysiological models will be available as
standalone models that can then be imported as appropriate into a new
composite model. Similarly for enzymes in metabolic pathways and
signalling complexes in signalling pathways, etc. Some examples of these
protein modules are:
Sodium/hydrogen exchanger 3 https://models.physiomeproject.org/e/236/
Thiazide-sensitive Na-Cl cotransporter
https://models.physiomeproject.org/e/231/
Sodium/glucose cotransporter 1
https://models.physiomeproject.org/e/232/
Sodium/glucose cotransporter 2
https://models.physiomeproject.org/e/233/
Note that in each case, as well as the CellML-encoded mathematical
model, links are provided (see Fig. 41) to the UniProt Knowledgebase
for that protein, and to the Foundational Model of Anatomy (FMA)
ontology (via the EMBLE-EBI Ontology Lookup Service) for information
about tissue regions relevant to the expression of that protein (e.g.
Proximal convoluted tubule, Apical plasma membrane; Epithelial cell
of proximal tubule; Proximal straight tubule). Similar facilities are
available for SMBL-encoded biochemical reaction models through the
Biomodels database [AYY].
Footnotes
Using PMR with OpenCOR
In addition to the PMR window for browsing public exposures directly in OpenCOR (PMR window) OpenCOR has the ability for users to directly create and access their workspaces in PMR.
Note
It is a feature of PMR that all data is persistent and permanent. As such, any workspaces created on the main instance of PMR (https://models.physiomeproject.org/) can not be deleted. For the purposes of teaching, we have an alternate instance of PMR (https://teaching.physiomeproject.org/) which is periodically cleared out and synschronised from the main instance. Using the teaching instance allows you to play around without the worry of things being permanent.
Register for a user account on the teaching instance of PMR.
In order to make use of the teaching instance of PMR, you must first have an account for that instance of the repository. For teaching purposes it is best to register a new account. This can be done by first opening this link in your browser: https://teaching.physiomeproject.org. Then click the Log in button on (shown in Fig. 42) then the registration form link.
After filling in the names and email fields and clicking Register you will receive an email inviting you to confirm and set a password. Once that is completed you can then log in. Clicking on My Workspaces will take you to a listing of all your workspaces and provides access to the the Workspace creation form.
Create a new workspace, in this example the title ‘Test workspace’ has been used.
The PMR Workspaces window
A window labelled PMR workspaces is available in OpenCOR (see Fig. 43). If it is not currently visible it can be selected via (or perhaps the Ctrl-space shortcut).
Set preferences.
Clicking the preferences button (Fig. 43) presents a Preferences dialog box with three settings: PMR instance, Name and Email. For the current purpose choose https://teaching.physiomeproject.org for the first and enter your name and email. These are used to identify you as the author of changes you submit back to the repository (view an example history).
Log into PMR from OpenCOR.
Before you can view private information or submit changes to PMR you must first log in to PMR from OpenCOR and grant OpenCOR permission to use your account. You accomplish this by clicking on the top right button in the PMR Workspaces window and then logging in with your new user name and password (created in step 1). Then grant access for OpenCOR to gain access to your PMR workspaces. The PMR workspaces window will then show all your workspaces, which should currently consist of the new workspace created in step 2. Note that using the same top right button you can log off - and when you next authenticate you will again be asked to grant access but this time without needing to login with your password.
Right clicking on the workspace name brings up a list of options for that workspace, the first being to view the workspace within PMR (in the web browser). Another option allows you to make a local copy of a workspace on your local disk - this will create a copy of the workspace on your local computer in which you are able to make changes.
Make a local copy of your test workspace
Using the option from the right-click menu on the workspace you created in step 2, clone the workspace to your PC. When doing this you will need to provide the folder in which you want to store the workspace contents - make sure you remember where this folder is!
Save a CellML model to your workspace.
A CellML file opened in OpenCOR (choose any model you have access to) can be saved () to the folder you created for the cloned workspace. Once you have saved a model you will see the file appear under the workspace’s folder in the PMR Workspaces window. Note that the file appears under the workspace with a red patch on the logo indicating the the file is not yet flagged to upload. To upload the file to PMR, you need to choose from the right-click menu on the workspace folder. This will ask you to provide a description of the change you would like to submit to PMR, and display all the differences you will be synchronising. When you now click the OK button, the changes will actually be submitted to PMR and you will see the file appear on the refreshed browser window. The file icon in the PMR Workspaces window will be shown without the red or green patch. Fig. 44 shows two CellML files that have been uploaded to PMR.
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 COMBINE 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. 33).
The recent versions of OpenCOR (since early 2016) support 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. 45.
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. 46 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
CellML provides a good technology to create, describe, and share definitions of mathematical models. SED-ML similarly provides a good technology to share reproducible descriptions of simulation experiments. Whenever possible, it is best to make use of these standard formats to ensure the models and simulations are Findable, Accessible, Interoperable, and Reusable.
Often in research projects, however, it is not always possible to describe the model and/or simulation that you need to perform in these declarative formats. It also doesn’t make sense to try and standardise extensions or modifications in such standards for potentially short-lived, one-off, research studies. Thus having access to a flexible scripting environment that works in concert with a standards-based tool like OpenCOR allows users to make use of standards when possible but with the flexibility to adapt as needed. OpenCOR supports this through the integration of a Python interpreter within the OpenCOR application.
Python-enabled versions of OpenCOR are now relatively mature, but still undergoing extensive user testing and implementation review. As such, this functionality is only available in special snapshot releases of OpenCOR available from: https://github.com/dbrnz/opencor/releases. In this part of the tutorial we are going to be using the 20 September 2019 snapshot of the Python-enabled OpenCOR. This particular release is distributed with the following Python packages and their dependencies: numpy, scipy, and matplotlib.
Python-enabled OpenCOR release can be installed as per the standard installation instructions. As this is an early release of the new functionality, it is best to use one of the compressed archive releases which you can extract locally rather than overwriting the stable system install. Once you have a Python-enabled release of OpenCOR your main OpenCOR window should look similar to that shown in Fig. 47.
OpenCOR application with default positioning of dockable windows including the Python console (right-side, middle). As described in Install and Launch OpenCOR the dockable windows can be rearranged as desired to suit your preferred layout.
In Python-enabled versions of OpenCOR the Python interpreter is embedded within the OpenCOR application. Which means that in order to access the OpenCOR functionality you must use that Python within the OpenCOR application rather than, for example, importing OpenCOR into your system Python. The Python console available in the OpenCOR graphical user interface handles this for you allowing a seamless user experience. However, often with Python-scripted simulation workflows it is nice to have the ability to run in a headless or batch mode. As such, Python-enabled versions of OpenCOR come with some command line scripts to help provide the user avoid the issues of making sure their Python scripts run using the correct Python interpreter.
In the top-level folder of your Python-enabled OpenCOR installation there is a script named run_python which will depend on your operating system - on Windows for example, it is called run_python.bat. Running this script without providing a Python script to execute will give you a standard Python console using the Python embedded inside the OpenCOR application:
C:\Users\andre\OpenCOR-2019-09-20-Windows> run_python.bat
Python 3.7.4 (default, Sep 20 2019, 18:29:34) [MSC v.1916 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>>
Providing a Python script will result in that script being interpreted by the interpreter embedded in the OpenCOR application:
C:\Users\andre\OpenCOR-2019-09-20-Windows> run_python.bat hello_world.py
Hello World!
Command line arguments can be provided as usual following the script to be executed.
Warning
Due to the use of a Python interpreter embedded in a graphical user interface, there can be some weirdness when trying to use UI toolkits from the command line, for example using matplotlib. This works within the OpenCOR graphical user interface, but will fail when running from the command line. Hence, it is best to currently use the command line version when working in a truly headless manner without the need for a graphical user interface.
There is another mode to make use of the Python-enabled versions of OpenCOR and that is to access this functionality via Jupyter notebooks. This is enabled via the run_jupyter helper script.
Todo
Write this section on Jupyter notebooks and OpenCOR.
As described above, the Python interpreter lives inside the OpenCOR application – making it difficult to access in order to install packages or modules that are not distributed with the Python-enabled versions of OpenCOR. To install packages using pip combined with the interactive Python console in the OpenCOR graphical user interface is the way to go here, as shown below.
Jupyter QtConsole 4.5.5
Python 3.7.4 (default, Sep 20 2019, 18:29:34) [MSC v.1916 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 7.8.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: !pip install [options] package
The scope and capabilities of the Python interface to OpenCOR is still being refined, but here we focus on use of the capabilities relevant to performing simulation experiments. Here we walk through the basic usage of using Python to interact with OpenCOR in performing a simulation experiment.
As with any Python script, we must first import the OpenCOR module to expose the functionality that we desire.
The main object that we are interested in dealing with is OpenCOR’s representation of a simulation. OpenCOR is able to generate a default simulation for a CellML model or to load a SED-ML document which defines the simulation experiment in detail. As the exposed simulation features are not yet complete, it is best to load a SED-ML document giving full control over the simulation settings. The output plots defined in the SED-ML will also be used when running the Python code via the interactive Python console in OpenCOR, but will be disregarded when running via the command line mode.
# for a local file
simulation = oc.openSimulation('path/to/cellml/or/sedml')
# OR for loading a remote file, e.g., from the model repository:
simulation = oc.openRemoteSimulation('URL/of/cellml/or/sedml')
# OR if using the OpenCOR GUI and models are already loaded
simulation = oc.simulation() # The model in the currently active tab
For a given simulation, the data object houses all the relevant information and pointers to the OpenCOR internal data representations.
And the data object allows us to define the interval of interest for this simulation experiment.
data.setStartingPoint(start)
data.setEndingPoint(end)
data.setPointInterval(pointInterval)
As in the OpenCOR graphical user interface, constant parameters and initial values for the state variables can also be set via the Python interface OpenCOR provides. When address specific variables in a model, they are mapped to Python dictionaries using key’s comprising of component_name/variable_name. This provides a method to uniquely identify all variables in a model.
# Set constant parameter values
data.constants()['key'] = value
# Set initial value for state variables
data.states()['key'] = value
Once you have the simulation defined that you would like to perform, it can be executed with the following.
If you are using the OpenCOR graphical user interface and have define plots for the current simulation experiment, then these will be displayed as usual during the execution of the simulation. The simulation results can also be used directly in the Python script as shown below.
# Access simulation results
results = simulation.results()
# grab a specific state variable results
r1 = results.states()['key'].values() # Numpy array
# grab a specific algebraic variable results
r2 = results.algebraic()['key'].values() # Numpy array
# access the full datastore representation of the simulation results
ds = results.dataStore()
# the dictionary or all result variables in the simulation
variables = ds.voiAndVariables()
# grab a the results for a given variable
r3 = variables['key'].values() # Python list of values
When continuing a simulation from an existing state, the default behaviour is to continue from the current state. The system can be reset to the initial state as shown below. As with using the OpenCOR graphical user interface, this includes resetting any parameters or initial values that you may have set via the GUI or the Python interface.
# Reset things if needed when re-running
simulation.resetParameters()
# clear any existing results
simulation.clearResults()
In this example, we use the simple ODE model introduced earlier in the tutorial. We will be using the Python console in the OpenCOR graphical user interface, working with the SED-ML loaded directly from the Physiome Model Repository. As we are using the OpenCOR application, you should see the user interface updating in response the to various Python commands. The following commands should be copy-pasted one at a time into the Python console to observe the behaviour.
import OpenCOR as oc
simulation = oc.openRemoteSimulation('https://models.physiomeproject.org/workspace/25d/rawfile/60ac9389285471a704f2f4be6e1a8ba5cbf45d1a/Firstorder.sedml')
data = simulation.data()
data.setStartingPoint(0)
data.setEndingPoint(10)
data.setPointInterval(0.1)
simulation.run()
# reset
simulation.resetParameters()
simulation.clearResults()
# change parameter values
data.constants()['main/b'] = 5
data.states()['main/y'] = 2
simulation.run()
# look at the simulation results
results = simulation.results()
y = results.states()['main/y'].values() # Numpy array
print(y)
ds = results.dataStore()
variables = ds.voiAndVariables()
y = variables['main/y'].values() # Python list of values
print(y)
a = variables['main/a'].values()
print(a)
In working through this example, you should be able to reproduce the results as seen in Fig. 7.
TensorFlow is a popular end-to-end open source machine learning platform in Python. Together with the Python-enabled OpenCOR capabilities and CellML itself, this opens up a new world of application of machine learning in computational physiology. This is a very new application that we are still actively developing, but here we give a brief demonstration that might help show what could be achieved.
The first step is to ensure that you have TensorFlow installed. As described above, Python packages need to be installed in the Python embedded inside OpenCOR. We are using here TensorFlow version 1.15, which can be installed using the OpenCOR Python console with the following command. (TensorFlow 2.0 will not work with this demonstration.)
In [1]: !pip install tensorflow==1.15
We have prepared a couple of Python scripts that you can use for this demonstration. The first is MPL.py, which is a TensorFlow-based script to construct a simple MLP (fully-connected feed-forward network or MultiLayer Perceptron) and trains it with a given dataset. The second is train-tf-model.py, which will first generate a set of training data using the O’Hara & Rudy cardiac electrophysiology model, which has been encoded in the CellML format as an extension of this model in the Physiome Model Repository. Both files should be downloaded into the same folder on your local machine.
Finally, in the OpenCOR Python console we need to make sure the plotting happens in-place rather than trying to bring windows. This is done by exectuing the following command in the OpenCOR Python console.
In [1]: %matplotlib inline
The train-tf-model.py script is the one that contains the definition of the workflow we are demonstrating here. It is easiest to open this file in your preferred Python editor and follow through the script, with the comments attempting to explain what is happening.
This script can be run in the OpenCOR Python console by first making sure the console is looking at the correct folder,
In [1]: %cd path/to/folder/with/downloaded/scripts
and then running the training script as follows.
In [1]: %run train-tf-model.py
All going well, this should result in something similar to Fig. 48.
You should now be able to play around with the training script to see what happens as you change, for example, the stimulation period or simulation duration.
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.
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
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.
