A local sensitivity analysis method for developing biological models with identifiable parameters: Application to cardiac ionic channel modelling

Future Generation Computer Systems 29 (2013) 591–598

Contents lists available at SciVerse ScienceDirect

Future Generation Computer Systems journal homepage: www.elsevier.com/locate/fgcs

A local sensitivity analysis method for developing biological models with identifiable parameters: Application to cardiac ionic channel modelling Anna A. Sher a,∗ , Ken Wang b , Andrew Wathen c , Philip John Maybank e , Gary R. Mirams a , David Abramson d , Denis Noble a , David J. Gavaghan e a

Department of Physiology, Anatomy and Genetics, University of Oxford, UK

b

Doctoral Training Centre, University of Oxford, UK

c

Department of Mathematics, University of Oxford, UK

d

Monash University, Australia

e

Oxford University Computing Laboratory, UK

article

info

Article history: Received 6 April 2011 Received in revised form 10 September 2011 Accepted 17 September 2011 Available online 15 October 2011 Keywords: Local sensitivity analysis Parameter identifiability Overparameterisation Cardiac modelling Ionic channels Singular value decomposition method

abstract Computational cardiac models provide important insights into the underlying mechanisms of heart function. Parameter estimation in these models is an ongoing challenge with many existing models being overparameterised. Sensitivity analysis presents a key tool for exploring the parameter identifiability. While existing methods provide insights into the significance of the parameters, they are unable to identify redundant parameters in an efficient manner. We present a new singular value decomposition based algorithm for determining parameter identifiability in cardiac models. Using this local sensitivity approach, we investigate the Ten Tusscher 2004 rapid inward rectifier potassium and the Mahajan 2008 rabbit L-type calcium currents in ventricular myocyte models. We identify non-significant and redundant parameters and improve the models by reducing them to minimum ones that are validated to have only identifiable parameters. The newly proposed approach provides a new method for model validation and evaluation of the predictive power of cardiac models. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Cardiac mathematical modelling is widely used to describe and predict physiological functions of the heart, as well as to reduce significantly the cost of experimental studies [1,2]. In mammalian ventricle, potassium and calcium currents are among the key ionic currents regulating the heart rhythm. The rapid delayed rectifier potassium current (IKr ) is important in the regulation of electrical activity of the heart, in particular the repolarization phase.1 Disturbances and mutations in IKr are associated with the prolongation of QT interval2 or with acquired Torsades de Pointes arrhythmias which are associated with sudden cardiac death.

∗

Corresponding author. Tel.: +44 0 1865 282 503. E-mail address: [email protected] (A.A. Sher).

1 The repolarization refers to the fall in the cardiac cell’s membrane potential that returns the membrane potential towards the resting potential. 2 The QT interval denotes the time interval between the start of the Q wave and the end of the T wave in the electrical cycle of the heart. 0167-739X/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.future.2011.09.006

Calcium ions play a significant role in the regulation of both electrical and mechanical activity of the heart [3]. One of the most important cardiac ionic currents is the L-type calcium current (ICaL ) as it enables the excitation–contraction coupling by allowing initial entry of Ca2+ inside the cell. Mutation or blockage of L-type calcium channels (LCCs) can lead to cardiovascular disorders (e.g. Timothy syndrome), causing cardiac arrhythmias and in some cases sudden cardiac death. Therefore, the use of mathematical ICaL and IKr models as a tool to explore the underlying mechanisms and cellular effects of ICaL and IKr can be extremely important, in particular when investigating drug induced arrhythmias. One of the common drawbacks of complex nonlinear mathematical cardiac models is overparameterisation. In this paper we present a new local sensitivity analysis method which allows to reliably assess parameter identifiability in a given model for a given set of experiments. Using the newly developed method, we (i) analyse the ICaL formulation in the Mahajan 2008 model [4], (ii) analyse the IKr formulation in the Ten Tusscher 2004 model [5], (iii) demonstrate the existing overparameterisation, and (iv) develop reduced models which can reproduce the key properties of IKr and ICaL currents, respectively.

592

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598

2. Methodology Local sensitivity analysis is an important tool for examining parameter identifiability3 in biological models [6,7]. Issues encountered with parameter estimation include model overparameterisation, and, specifically, the following challenges are faced: (i) varying parameter values locally produces non-significant shifts in the model output, (ii) parameters may be linearly dependent, so that removing one redundant parameter and at the same time varying the other linearly dependent parameters accordingly yields the same model output. 2.1. Parameter identifiability and local sensitivity analysis To investigate the significance and redundancy of model parameters, local sensitivity methods can be used. The local effect of a parameter Pj on the output variable yi at time t can be described by sensitivity factor Jij Jij (t ) =

∂ yi (t ) ∂ Pj

(1)

with sensitivity matrix J being constructed such that each column describes the sensitivity of all outputs with respect to a particular parameter. Existing approaches include, for example, Eigenvalue Decomposition (EVD) and QR based4 algorithm developed by Fink et al. [6], where (1) the EVD is performed on J T J (or on J if J is a square matrix) and eigenvalues normalised with respect to the maximum eigenvalue are recorded, and (2) the QR decomposition is performed on J T J (or on J if J is a square matrix) to obtain a permutation vector E, so that the final score is calculated by associating the normalised eigenvalues with the parameters according to E. Such existing methods may provide partial insights into the significance of the parameters. However, even though this approach helps in categorising parameters according to whether or not these are identifiable or non-identifiable, it cannot determine the cause of non-identifiability. In particular, such approach is unable to identify redundant parameters in an efficient manner. To identify redundant and non-significant parameters in an efficient manner, we propose a new Singular Value Decomposition (SVD) [8] based approach. Newly developed SVD based parameter identifiability algorithm: While the above EVD approach focuses on identifying the most significant parameters, the newly proposed SVD approach focuses on identifying the redundant and least significant parameters. Specifically, information about the redundant parameters can be determined by investigating linear dependence of the columns of J. In particular, in the proposed method, (a) Singular value decomposition is performed on J such that J = USV ∗ ,

(2)

where U is an m by n unitary matrix, S is a diagonal matrix, V is an n by n unitary matrix and V ∗ denotes the conjugate transpose of V ; (b) Linearly dependent parameter subset(s) are identified by investigating vectors vi which span the null space of J, so that J ∗ vi = 0,

(3)

where vi defines the subset of the singular vectors that correspond to vanishing values of J.

3 A parameter is defined as identifiable if it is both significant and non-redundant. 4 The QR refers to the QR decomposition of a real matrix in which we want to compute the eigenvalues; by definition, Q denotes an orthogonal matrix and R represents an upper triangular matrix.

Fig. 1. The Ten Tusscher 2004 IKr formulation with highlighted linearly dependent parameters that are identified by the newly developed SVD method.

The final score for each parameter is calculated by summing the contributions of the parameter in each of the parameters’ directions in the space spanned by the sensitivity matrix, represented by entries in V and associated with corresponding singular values. 2.2. Cardiac modelling and numerical simulations For parameter identifiability studies we use (1) the Ten Tusscher 2004 IKr Hodgkin Huxley type model of human ventricular myocyte rapid inward rectifier potassium current [5], and (2) the Mahajan 2008 ICaL 7-state Markov Model of rabbit ventricular myocyte L-type Calcium Current [4]. We choose the Ten Tusscher 2004 human IKr [5] model as it is amongst the most widely used IKr models and is also one of the models analysed by [6] using the EVD method. We choose the Mahajan 2008 ICaL model as it is one of the most advanced ICaL models, being developed using both training and validation data sets. Specifically, it was derived based on experimental data when fitting the transition rates as well as other parameters related to voltage- and calcium-dependent inactivation of ICaL . Table 2 summarises the training and validation protocols. All models and methods are coded using MATLAB R2008a. The ordinary differential equations are solved using ode15s, with the maximum relative error tolerance of integration set to 10−6 . To examine parameter identifiability (and next, if required, to reduce the model to include only identifiable parameters), we first need to identify model parameters of interest as well as experimental protocols (both training and validation sets) used to estimate these parameters (see section below). 2.3. Experimental and simulation protocols and choice of model parameters In this work, we examined the influence of 14 and 30 model parameters on the shape of IKr and ICaL currents respectively using local sensitivity analysis methods: EVD and SVD. Protocols and parameters investigated are summarised below and in Figs. 1 and 2 (upper panel) and Tables 1 and 2. The comparison between two different local sensitivity methods is performed to assess their respective strengths, and the results are summarised in Section 3. In the case of IKr we perform sensitivity studies using the EVD and SVD methods to assess 14 parameters of the Ten Tusscher 2004

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598

60 40 20 V(mV)

0 -20 40 -60 -80 100 0

5

10

15

10

15

Time (s) 50 Orginal parameters Manually refitted parameters 40

IKr(A/F)

30

20

10

0

–10

0

5

593

from a holding potential of +60 mV (for 1000 ms) to −100 and −80 mV for 10 ms each followed by steps to +40 mV for 100 ms, to −60 mV for 200 ms and back to the holding potential. In the case of ICaL we assess the parameter estimation of voltagedependent inactivation (VDI) and calcium-dependent inactivation (CDI) of ICaL related parameters. We use the training experimental protocols as described in [4]. In particular, we consider voltage clamp protocols where standard rectangular pulse protocols are delivered as step pulses from 30 to +30 mV in 10 mV increments from a holding potential of −80 mV, holding for 300 ms at each step, under various conditions, namely (i) control, (ii) depleted sarcoplasmic reticulum (SR), (iii) 1.8 mM Ba, where Barium (Ba2+ ) is introduced instead of extracellular calcium of 1.8 mM Cae . As in [4], we use separate protocols to estimate CDI and VDI transition rates. Specifically, to fit voltage-dependent inactivation ICaL transition rates, we fit VDI related parameters to reproduce the ICaL recorded under voltage clamp with SR-depleted and 1.8 mM Ba. To fit CDI related parameters, both control voltage clamp and SR-depleted (with 1.8 mM Ca2+ ) protocols are used. To mimic the estimation procedure in [4], we split parameters of interest into the following three subsets: (a) parameters involved only in the CDI pathway (pure CDI parameters), (b) parameters involved only in the VDI pathway (pure VDI parameters), and (c) parameters involved in both pathways (shared CDI and VDI parameters). The list of parameters used for estimation5 is summarised in Table 1. ICaL related parameters which can be uniquely determined by the independent fitting from the above set of experiments are excluded from the list of parameters of interest shown in Table 1. In particular, Pca (LCC open probability), gca (maximum conductance of LCC), as well as R(V ) related parameters are not used for sensitivity studies in this work, and hence not shown in Table 1, as these are estimated using an independent set of experiments such as the density of LCCs, maximum conductance rate of LCC, and voltage clamp recovery protocol, respectively.

Time (s) Fig. 2. UPPER PANEL: The Fink and Noble [6] short voltage clamp protocol. LOWER PANEL: The Ten Tusscher 2004 IKr current under the Fink and Noble [6] short voltage clamp under control conditions (blue) and with redundant parameters removed (red), i.e. where p(8) and p(14) are removed by refitting the original parameters p(5) = 450, p(8) = 6, p(11) = 2 and p(14) = 1.12 to become p(5) = 2700, p(8) = 1, p(11) = 3.36 and p(14) = 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

IKr model that describe voltage dependent activation, inactivation, deactivation and reactivation of the IKr channel [5] (Fig. 1). Parameters such as the Nernst potential for potassium (EK ), maximum conductance of the IKr channel (gKr ) and extracellular potassium concentration (Ko ) are not included in the set of 14 parameters analysed since these can be estimated using an independent set of experiments such as the concentration of potassium ions inside and outside the cell, density of the IKr channels and their maximum conductance rate. We investigate the sensitivity of IKr to voltage dependent gating parameters under application of the Fink and Noble [6] voltage clamp protocol (Fig. 2, upper panel). The latter is derived to include the information content of a normal range of voltage-clamp steps necessary to characterise voltage-gated ion channel dynamics and is defined as follows [6]: (i) voltage-dependent activation and inactivation stage: step pulses from a holding potential of −80 mV (for 600 ms) to +60, +40 and +20 mV for 1000 ms, and back to −80 mV after each pulse; (ii) voltage-dependent deactivation stage: step pulses from a holding potential of +60 mV (for 1000 ms) to −100, −70 and −40 mV for 1000 ms and back to the holding potential after each pulse; (iii) voltage-dependent reactivation stage: step pulses

3. Results and discussion In our sensitivity studies outlined below, we use protocols summarised in Section 2.3 to perform the following analysis. First, we perform sensitivity analysis using the training protocols to determine whether or not the IKr voltage related gating parameters and the ICaL CDI, VDI and shared parameters are identifiable, and compare results of applying EVD and SVD methods, highlighting the advantages and disadvantages of each. The sensitivity results are verified to check if the non-identifiable parameters can be removed and whether or not the reparameterised newly proposed reduced IKr and ICaL models reproduce the original respective models’ behaviour under application of training protocols. Finally, we validate the newly proposed reduced IKr and ICaL models to contain only identifiable parameters and, in the case of ICaL , we also verify whether the reduced ICaL model is able to reproduce the results of the original ICaL model under application of validation protocols. Tables 3 and 4 show the sensitivity results of EVD and SVD methods for IKr and ICaL models respectively. Importantly, the newly developed SVD method provides information both on

5 Often parameters presented as model parameters in tables inside the papers do not form the full list of parameters that the model contains and which were estimated based on the experimental data. For instance, as shown in Fig. 4 in equation for TCa , the value of 0.1 is not listed as a parameter in [4], yet it is one of the parameters which was fitted by [4] using the training data sets. Thus, one has to carefully define the list of model parameters for estimation.

594

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598 Table 1 ICaL related parameters in the Mahajan 2008 model. Category

ICaL related parameters

Only VDI related

k2t s1t k1t const-tauba

Only CDI related

cpt cat k2 tca const-Tca2 const-S const-k1 power-Tca const-Tca1 power-fca const-fca

VDI and CDI shared

taupo r1 r2 s6 vx sx tau-3 vyr const-pr const-ps const-poinf syr vy sy const-k3

non-significant parameters and linearly dependent parameter subset(s) (hence, redundant parameter(s)). As discussed in Section 2.1, the EVD method can group parameters as identifiable and nonidentifiable. Yet, the EVD method, as opposed to the SVD method, cannot distinguish the cause of the non-identifiability of parameters, i.e. whether or not the non-identifiable parameters are nonsignificant or linearly dependent. 3.1. IKr sensitivity studies 3.1.1. IKr model and parameter identifiability Table 3 shows the sensitivity results of EVD and SVD methods for the IKr model. Two sets of linearly dependent parameters, (1) p5 and p8 and (2) p11 and p14 , are identified by the SVD method. The EVD method, however, only finds two non-identifiable parameters, p8 and p11 . This result shows that not only the EVD method cannot determine the cause of the parameter non-identifiability, i.e. whether parameters are non-significant or linearly dependent, but also the method cannot identify the full set of the linearly dependent parameters. Thus, the EVD method provides limited information on the non-identifiable parameters.

3.1.2. Reduced IKr model The mathematical formulation of the IKr Ten Tusscher 2004, Hodgkin–Huxley type, model is presented in Fig. 1 where the nonidentifiable parameters obtained by the SVD method (Table 3) are highlighted by arrows. To improve the IKr formulation, we remove redundant parameters identified by the SVD method. Specifically, we exclude parameters p(8) and p(14) and manually refit the original IKr model (Fig. 2). The original parameters p(5) = 450, p(8) = 6, p(11) = 2 and p(14) = 1.12 are refitted to become p(5) = 2700, p(8) = 1, p(11) = 3.36 and p(14) = 1. The reduced IKr model contains 12 VDI dependent parameters instead of 14 parameters of the original model. Such reduction in the number of parameters would yield a significant reduction in computational costs, for instance, when this IKr ionic model is used for tissue cardiac simulations. Fig. 2 shows the Ten Tusscher 2004 IKr current under the Fink and Noble [6] short voltage clamp under control conditions (blue) and with redundant parameters removed (red). As expected, the reparameterised reduced IKr model shows identical output to the control case. 3.1.3. Reduced IKr model and parameter identifiability Applying SVD analysis to the reduced IKr model, we validated that the new IKr model contains no non-identifiable parameters. Applying EVD analysis to the reduced IKr model also yielded no non-identifiable parameters. This verification by the SVD method (to check whether there are any non-significant or linearly dependent parameters in the new IKr model) shows that the reduced IKr model contains a minimum set of parameters required to describe the whole-cell rapid delayed rectifier potassium current. 3.2. ICaL sensitivity studies 3.2.1. ICaL model and parameter identifiability Table 4 shows the sensitivity results of EVD and SVD methods for the ICaL model. Further detailed comparison of the sensitivity results obtained by the EVD and SVD methods is presented in Table 5, where the results are shown for four parameters (Tca, constTca2, constTca1, cpt) which were identified as non-significant or linearly dependent by the newly developed SVD method. As discussed in Section 2.1, while the EVD method can group parameters as identifiable and non-identifiable, it cannot distinguish the cause of the non-identifiability of parameters, i.e. whether or not the non-identifiable parameters are nonsignificant or linearly dependent. Moreover, as shown in Table 4, the EVD method identified a larger number of non-identifiable parameters than the SVD method. However, these results are misleading, since the normalisation according to the largest eigenvalue (inherent feature of the EVD method, Section 2.1) may lead to false results sometimes as the parameter score may not accurately correlate with the parameters’ significance. For example, results in Fig. 3 illustrate that some of the parameters

Table 2 Experimental protocol and corresponding simulated protocol used for the Mahajan 2008 ICaL sensitivity analysis. Experimental protocols in [4]

Simulations performed

Training protocols: Voltage step clamp under normal conditions Voltage step clamp under SR depleted conditions where SR uptake and release are blocked by drugs Voltage step clamp under i.e. Barium protocol: Barium influx instead of Calcium to remove CDI of ICaL

Original model parameters are used Flux from SR is zeroed; Submembrane Ca2+ is clamped to its initial value CDI transition rates are turned to zero

Validation protocols: AP clamp

Original model parameters are used

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598

595

Table 3 Parameter identifiability results of the EVD and the newly developed SVD algorithms for the IKr current under the Fink and Noble [6] short voltage clamp protocol. Parameter category

VDI related

Non-identifiable parameters found by each method: Existing Eigenvalue Decomposition algorithm

Newly developed SVD algorithm

p8 p11

Linearly dependent:p5 and p8 , p11 and p14 Non-significant: N/A

Table 4 Parameter identifiability results of the EVD and the newly developed SVD algorithms for the ICaL current under training set protocols from Table 2. Parameter category

Non-identifiable parameters found by each method: Existing Eigenvalue Decomposition algorithm

Newly developed SVD algorithm

Only CDI related

power-fca power-Tca k2 tca constTca1 const-k1 constTca2 const-s1

Linearly dependent: cpt tca constTca1 Non-significant: constTca2

Only VDI related

None

None

VDI and CDI shared

r2 syr sy tau3 sx

None

Table 5 Local sensitivity results of the EVD and SVD based algorithms for the ICaL Mahajan 2008 rabbit model for the CDI related parameters identified as non-significant or linearly dependent by the SVD method in Table 4. Parameter

Existing Eigenvalue Decomposition algorithm

Newly developed SVD algorithm

constTca2

Low score

Non-significant Non-redundant

Tca

Medium/High score

Highly significant Linearly dependent

constTca1

Low score

Medium significant Linearly dependent

cpt

High score

Medium/low significant Linearly dependent

found as non-identifiable by the EVD method are actually significant. Specifically, as demonstrated in Fig. 3, parameters power-Tca and power-fca, which were found to be non-identifiable by the EVD method, actually have a significant effect on the shape of ICaL , hence the EVD method is misleading in suggesting that these parameters are non-identifiable. 3.2.2. Reduced ICaL model with identifiable parameters The mathematical formulation of the ICaL 7-state Markov model is shown in Fig. 4, where the non-identifiable parameters obtained by the SVD method (Table 5) are highlighted by arrows and boxes. To improve the ICaL model, we reduce the model to exclude the non-identifiable parameters. Specifically, the reduced ICaL model is obtained by removing non-significant and redundant parameters, where, if any redundant parameters were found, the reduced model would be a reparameterised version of the original model.6

6 When removing redundant parameter(s) from the original model, the remaining linearly dependent parameter(s) are refitted accordingly in the newly obtained reduced model to accommodate for the removal of the redundant parameter(s).

As demonstrated in the upper panel of Fig. 5, removal of parameter constTca2, which was identified as non-significant by the SVD method, indeed yields no effect on the ICaL current under the voltage clamp experiment. Moreover, similarly no effect is seen on the ICaL shape when constTca2 is removed under the application of voltage clamps under SR depleted conditions and (or) Barium protocols described in Table 2 (data not shown). Fig. 5 shows results of the removal of the redundant parameter cpt, which is chosen to be the one to be removed as among the three linearly dependent parameters it is the least significant (see Table 5). Removal of cpt parameter (red curve, lower panel of Fig. 5) without reparametrization of the model affects the ICaL shape, which is consistent with the SVD method results, as the parameter was identified as linearly dependent but has not been found to be non-significant. As expected, however, after reparameterising the model upon the removal of the redundant cpt parameter, the resulting model produces ICaL results identical to the control case (green curve, lower panel of Fig. 5) upon application of voltage clamp protocols which were used to identify the parameter identifiability. Similar results are obtained when applying voltage clamps under SR depleted conditions and (or) Barium protocols described in Table 2 (data not shown).

596

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598

0

parameters, which also indicates that the model of a given formulation contains a minimum set of parameters required to describe the ICaL recordings under the given voltage clamp protocols summarised in Table 2.

ICaL (A/F)

-0.05

–0.1

–0.15 control varied parameter power-Tca

–0.2

0

500

1000 Time (ms)

1500

2000

0

ICaL (A/F)

-0.05

–0.1

3.2.3. Validation of the reduced ICaL model Previous sections focused on the training experimental protocols in order to identify whether or not the model is overparameterised, and to reduce the model to a minimum one that contains only identifiable parameters when tested against the training data. In this section, we validate the new reduced ICaL model using validation protocol sets to confirm that the newly proposed model is able to reproduce the key properties. Also, the success of testing the validation protocol sets would also be indicative of whether or not the choice of training protocols is reasonable. As demonstrated in Fig. 6, applying the AP clamp, which is the validation protocol used by the Mahajan 2008 [4] (see Table 2), yields same results whether the original (blue curve) or the new reduced (red curve) ICaL model is used. Further, the model is successfully validated for the AP clamp protocol applied under (i) SR depleted and (ii) Barium protocol conditions (data not shown). These validation studies show that the new reduced model can be used instead of the original ICaL model. 4. Conclusion

–0.15 control varied parameter power-fca

–0.2

0

500

1000 Time (ms)

1500

2000

Fig. 3. ICaL plot in normal cell voltage clamp before and after varying power-Tca (4 to 1) (top plot), and power-fca (3 to 1) (bottom plot).

Next, applying SVD analysis to the reduced ICaL model, we validate that the new model does not contain non-identifiable

In this paper we present a new parameter identifiability method which is more powerful than the previously developed approaches. Applying the newly developed method to the Ten Tusscher 2004 IKr and Mahajan 2008 ICaL models allowed us to obtain minimal respective models which we validated to contain only identifiable parameters. The results presented are also important as they illustrate that even in the case of models constructed using advanced parameter estimation approaches (e.g. the Mahajan 2008 ICaL model developed based on training and validation protocols [4]) there is still a possibility of having an overparameterised system.

Fig. 4. The Mahajan 2008 ICal 7-state Markov model formulation with highlighted non-significant and linearly dependent parameters that are identified by the newly developed SVD method.

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598

0

The limitation of the developed method is that at present the reparameterisation step is done manually. In future we will automate this procedure; specifically, the newly developed method will be extended to improve the refitting procedure, so that nonlinear minimization methods are used to identify a new set of parameter values in reparameterised models.

ICaL (A/F)

-0.05

-0.1

References -0.15 control non-significant parameter (constTCa2) removed

-0.2

0

500

1000

1500

2000

Time (ms) 0 -0.1 -0.2 ICaL (A/F)

597

-0.3 -0.4 -0.5 control

-0.6 -0.7

[1] D. Noble, A. Varghese, P. Kohl, P. Noble, Improved guinea-pig ventricular cell model incorporating a diadic space, ikr and iks, and length- and tensiondependent processes, The Canadian Journal of Cardiology 14 (1) (1998) 123–134. [2] T.R. Shannon, F. Wang, J. Puglisi, C. Weber, D.M. Bers, A mathematical treatment of integrated Ca dynamics within the ventricular myocyte, Biophysical journal 87 (5) (2004) 3351–3371. [3] D.M. Bers, Excitation-Contraction Coupling and Cardiac Contraction Force, Kluwer Academic, 2001. [4] A. Mahajan, Y. Shiferaw, D. Sato, A. Baher, R. Olcese, L.-H. Xie, M.-J. Yang, P.-S. Chen, J.G. Restrepo, A. Karma, A. Garfinkel, Z. Qu, J.N. Weiss, A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates, Biophysical Journal (2008). [5] K.H.W.J.ten Tusscher, D. Noble, P.J. Noble, A.V. Panfilov, A model for human ventricular tissue, American Journal of Physiology 286 (2004) H1573–H1589. [6] M. Fink, D. Noble, Markov models for ion channels: versatility versus identifiability and speed, Physical and Engineering Sciences 367 (1896) (2009) 2161–2179. [7] D.M. Hamby, A review of techniques for parameter sensitivity analysis of environmental models, Environmental Monitoring and Assessment 32 (2) (1994) 135–154. [8] G.H. Golub, C.F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, 1996.

refitted removed cpt parameter without refitting

0

200

400

600

800

1000

Time (ms) Fig. 5. UPPER PANEL: The Mahajan 2008 ICal current under control conditions and with non-significant parameter constTca2 removed. The protocol is ICal recorded under voltage clamp from −80 mV to −20, −10, 0, 10, 20, 30 mV and depleted SR conditions. LOWER PANEL: The Mahajan 2008 ICal current under control conditions and with redundant parameter cpt removed. The protocol is ICaL recorded under voltage clamp from −80 mV to −20, −10, 0, 10, 20, 30 mV and depleted SR conditions.

0 -0.05

Anna A. Sher received her B.Sc. degree in Physiology and Mathematics from McGill University, Canada in 2002. From there she moved to the University of Oxford where she continued to pursue her interests in mathematical biology. At the University of Oxford she gained an M.Sc. in Applied and Computational Mathematics in 2003, and a D.Phil. in Computational Biology in 2008. During her M.Sc. Anna worked on electrochemical problems related to biosensors, modelling Frequency Domain AC Voltammetry. During her D.Phil. she investigated the role of the Na+ –Ca2+ exchanger in the regulation of calcium dynamics in cardiac cells during heart failure using mathematical modelling and numerical simulations. Currently she works as a postdoctoral researcher at the Cardiac Electrophysiology Group at the Department of Physiology, Anatomy and Genetics at the University of Oxford modelling drug action in ventricular myocytes at the cellular and sub-cellular levels to assess drug-induced pro-arrhythmic risk.

ICaL (A/F)

-0.1 -0.15

Ken Wang is a D.Phil. student in the Computing Laboratory at Oxford University. She holds a BA degree in Mathematics from the University of Oxford. She works on action potential heterogeneity in rabbit left ventricle using both experiments and mathematical modelling.

-0.2 -0.25 control reduced model

-0.3 -0.35

0

100

200

300

400

500

Time (ms) Fig. 6. The AP clamp protocol applied to the reduced Mahajan model where the AP is generated from the Mahajan original model. The ICaL from the original Mahajan model (blue) and the reduced model (red) are plotted here. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Thus, the method presented here is of importance when choosing the models to work with as well as when developing the models, as it allows to verify which parameters are non-significant and/or linearly dependent.

Andrew Wathen is Reader in Numerical Analysis at Oxford University. His research is in the general area of Numerical Linear Algebra and Partial Differential Equations, in particular on preconditioning and iterative methods which are of particular relevance for large scale scientific computing applications. His work on block preconditioners for saddle-point systems has impacted a number of application areas: it has led for example to fast solvers for incompressible flow simulation and a significant class of methods of general utility in Optimization which emanate from his introduction of ‘constraint preconditioning’.

598

A.A. Sher et al. / Future Generation Computer Systems 29 (2013) 591–598 Philip John Maybank is a D.Phil. student in the Oxford University Computing Laboratory at Oxford University. He holds a B.Sc. in Mathematical Sciences from Durham University, and an M.Phil. in Computational Biology from Cambridge University. He works on algorithms for model reduction using mathematical techniques from functional analysis.

Gary R. Mirams received his M.Math in mathematics with engineering from the University of Nottingham (2004). He went on to study at the University of Oxford’s Life Science Interface Doctoral Training Centre (2004–2005), before returning to Nottingham where he gained a Ph.D. in Mathematical Biology (2008) with a thesis entitled ‘‘Subcellular Phenomena in Colorectal Cancer’’. He is interested in automating model development, and is involved in the development of the Chaste computational biology software (http://www.comlab.ox.ac.uk/chaste/). He is currently working on the computational prediction of drug side-effects on the electrical activity of the heart, in conjunction with pharmaceutical company partners. David Abramson has been involved in computer architecture and high performance computing research since 1979. Previous to joining Monash University in 1997, he has held appointments at Griffith University, CSIRO, and RMIT. At CSIRO he was the program leader of the Division of Information Technology High Performance Computing Program, and was also an adjunct Associate Professor at RMIT in Melbourne. He served as a program manager and chief investigator in the Co-operative Research Centre for Intelligent Decisions Systems and the Co-operative Research Centre for Enterprise Distributed Systems. Abramson is currently an ARC Professorial Fellow; Professor of Computer

Science in the Faculty of Information Technology at Monash University, Australia, and Science Director of the Monash e-Research Centre.

Denis Noble CBE, FRS, FMedSci, is Emeritus Professor of Cardiovascular Physiology. He now directs the computational physiology research group. He was the first to model cardiac cells (in two papers in Nature in 1960) and has published over 350 research papers. He is one of the leaders of Systems Biology and has written the first popular book on Systems Biology, The MUSIC of LIFE (OUP, 2006).

David J. Gavaghan received his undergraduate degree in Mathematics from Durham University in 1986. From there he moved to the University of Oxford where he gained an M.Sc. in Numerical Analysis and Mathematical Modelling in 1987 and a D.Phil. in the Development of Parallel Numerical Algorithms in 1991. Since then he has been working in the field of Mathematical and Computational Modelling and has established and heads the Computational Biology Group which is based principally within the Department of Computer Science. In 2002 he was awarded substantial funding from the EPSRC to establish the Life Science Interface Doctoral Training Programme. Further programmes in Systems Biology (funded in 2007) and an industrial programme in Systems Approaches to Biomedical Science (funded in 2009) have since been added, and the Doctoral Training Centre is now based in the Rex Richards Building in the heart of the University’s South Parks Road Science Campus. In October 2004 he was appointed Professor of Computational Biology within the Computing Laboratory (now Department of Computer Science). His research interests are in the mathematical modelling of physiological and biological systems, and in the development of robust approaches to the development of computational science software that will result in fully tested and reliable (open source) codes.