Parameter estimation in cardiac ionic models

ARTICLE IN PRESS

Progress in Biophysics & Molecular Biology 85 (2004) 407–431

Parameter estimation in cardiac ionic models Socrates Dokos*, Nigel H. Lovell Graduate School of Biomedical Engineering, University of New South Wales, Sydney 2052, NSW, Australia

Abstract We examine the problem of parameter estimation in mathematical models of excitable cell cardiac electrical activity using the well-known Beeler–Reuter (1977) ionic equations for the ventricular action potential. The estimation problem can be regarded as equivalent to the accurate reconstruction of ionic current kinetics and amplitudes in an excitable cell model, given only action potential experimental data. We show that in the Beeler–Reuter case, all ionic currents may be reasonably reconstructed using an experimental design consisting of action potential recordings perturbed by pseudo-random injection currents. The Beeler–Reuter model was parameterised into 63 parameters completely defining all membrane current amplitudes and kinetics. Total membrane current was fitted to model-generated experimental data using a ‘data-clamp’ protocol. The experimental data consisted of a default action-potential waveform and an optional series of perturbed waveforms generated by current injections. Local parameter identifiability was ascertained from the reciprocal condition value (1/l) of the Hessian at the known solution. When fitting to a single action potential waveform, the model was found to be over-determined, having a 1/l value of B3.6e
14. This value improved slightly to B1.4e
10 when an additional 2 perturbed waveforms were included in the fitting process, suggesting that the additional data did not overly improve the identifiability problem. The additional data, however, did allow the accurate reconstruction of all ionic currents. This indicates that by appropriate experimental design, it may be possible to infer the properties of underlying membrane currents from observation of transmembrane potential waveforms perturbed by pseudo-random currents. r 2004 Elsevier Ltd. All rights reserved. Keywords: Parameter estimation; Cardiac ionic models; Non-linear fitting

*Corresponding author. Tel.: +61-293853929; fax: +61-296632108. E-mail address: [email protected] (S. Dokos). 0079-6107/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.pbiomolbio.2004.02.002

ARTICLE IN PRESS 408

S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

1. Introduction Over several decades, mathematical models of cardiac excitable cell electrical activity have provided important integrative insights into the nature of ionic currents underlying various phases of the cardiac action potential. These models have closely followed the framework pioneered by Hodgkin and Huxley (1952) in their model of a squid axon, consisting of a non-linear system of ordinary differential equations (ODEs). In this framework, whole-cell membrane or individual channel currents may be experimentally characterised using a voltage-clamp protocol, with individual current components teased apart using pharmacological intervention. Such experiments are able to yield quantitative characterisations of individual ionic current kinetics and conductances, such that this data may, in principle, be readily incorporated into mathematical models for the computational reconstruction of action potential activity. In practice, however, the quantitative data provided by these experiments is subject to a high degree of biological variability from preparation to preparation. This may be due, for example, to inherent differences in ionic channel expression among cells. Cardiac cells are known to exhibit spatial heterogeneity in action potential waveshape in many regions of the heart (Liu et al., 1993; Feng et al., 1998). From an experimental view, there are significant technical issues associated with the recording of small currents from beating cells, applied pharmacological agents used may not be sufficiently specific or selective for the purposes of computational modeling, and the absence of specific experimental data on an ionic current often necessitates the amalgamation of data from different species into a model. All of the above suggest that models assembled from incomplete and/or inaccurate data may result in limited predictive power (Noble and Bett, 1993; Hunter et al., 2001]. It is therefore highly desirable, arguably even necessary, to experimentally validate a given cardiac ionic model against data obtained from a specific cardiac cell. For cardiac ionic models, a particular problem with such cell-specific validation is the large number of parameters inherently present. Large numbers of parameters are likely to result in over-determined systems for which unique identification is not possible. In this article, we examine the problem of parameter identifiability in a well-known model of cardiac ventricular electrical activity, the Beeler and Reuter (1977) (BR) model, consisting of 63 parameters. Plots of the action potential, underlying membrane currents and cytosolic calcium concentration of this model are shown in Fig. 1. A full description of the BR equations, along with the specific parameterisation adopted in this study, is given in the appendix. Ionic models of cardiac action potentials electrically model the cell membrane as an equivalent capacitance in parallel with several transmembrane conductance paths, each describing an ionic current or exchanger. In space-clamped mode, the total ionic current (itot) is equal and opposite to the capacitive current, so that the transmembrane potential (Vm) is governed by the first-order differential equation dVm itot ¼
dt Cm

ð1Þ

where Cm is the membrane capacitance. The individual ionic currents, which sum to form itot, are often dynamic functions of Vm itself. They may be described by additional ODEs which characterise their kinetic properties. The resulting system of ODEs leads to a complex high-order system which must be solved by numerical integration. For the experimenter, the electro-

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

409

40 20

E

mV

0 -20 -40 -60 -80 -100 100 i stim

µA/cm

2

0 -100

iNa

-200

2 µA/cm

5

iK1 ix1

0

is -5 6 [Ca] i

µM

4 2 0

0

50

100

150

200

250 ms

300

350

400

450

500

Fig. 1. Standard BR model using default values for all parameters. The panels illustrate the transmembrane action potential (top), underlying membrane and external stimulus currents (middle panels) and the intracellular calcium concentration (lower panel). The stimulus current was set at 40 mA/cm2, duration 1 ms. In a typical experiment, the experimenter has access only to the Vm data in the top panel. The identification problem is to determine the conductance characteristics and kinetic properties of the hidden ionic currents, given that only Vm can be directly observed.

physiological problem of interest is to infer the individual nature and properties of the underlying ionic currents which make up itot in (1). Mathematically, this is equivalent to estimating model parameters which define current kinetics and/or other sub-cellular processes. In real cardiac cells, the experimenter can observe only Vm directly, and must therefore infer the behaviour of the

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

410

underlying ionic currents by indirect means, usually by injecting an external feedback current to hold Vm at a preset constant level (the ‘‘voltage-clamp’’ method). In the present article, we explore the alternative possibility of estimating model parameters using a least-squares fit approach to experimental data generated in the absence and/or presence of external pseudo-random perturbation currents. The advantage of this approach is the that it opens up the possibility of formulating, for the first time, cell-specific reconstructions of underlying ionic mechanisms using only a relatively simple set of experimental data generated from current injections.

2. Methods 2.1. Choice of parameterisation The equations defining the complete BR model of the ventricular action potential, simulated in Fig. 1, were parameterised into a system of 8 first-order ODEs employing 63 adjustable parameters. A complete description of model equations and adopted parameterisation is provided in the appendix. Where possible, all constants and coefficients appearing in any equation of the original BR model were selected as parameters to be optimised. These parameters included all scaling factors specifying membrane current amplitudes and fully activated current–voltage characteristics, as well as all coefficients of the individual rate equations. In most cases, these rate equations were rewritten into an equivalent form, in order to remove any obvious dependencies amongst parameters. As an example of the latter, the rate equation specifying the opening rate of the x1 gating variable (ax1 ) is given in the original BR model as ax1 ¼

C1 eC2 ðVm þC3 Þ ; eC4 ðVm þC3 Þ þ C5

ð2Þ

where C1 yC5 are specific parameters in the original description. However, the five parameters in (2) are not independent, since the equation can also be expressed in the following equivalent form employing only 4 parameters: ax1 ¼

p1 ep2 Vm ; ep3 Vm þ p4

ð3Þ

where p1  C1 eC3 ðC2
C4 Þ ; p2  C2 ; p3  C4 and p4  C5 e
C3 C4 : Hence, (3) instead of (2) was used in this study. Similar modifications were also applied to other rate formulations in the BR equations. Another source of parameter dependency lies with Eq. (1). Here, the total membrane current is scaled by a factor of Cm. As a result, any other parameters which set the scale of individual membrane currents will not be independent of Cm. We address this problem by fixing Cm at its default value, allowing the other current amplitude parameters to freely vary. This is equivalent to simply specifying the magnitude of all ionic currents per unit of cell membrane capacitance.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

411

2.2. Data-clamp technique In order to experimentally estimate the values of ionic model parameters, we assume that the experimenter has access only to the transmembrane action potential waveform (e.g. upper panel in Fig. 1). Because of this restriction, it is necessary to numerically integrate the full set of dynamic model equations for each candidate parameter set, since the behaviour of the underlying ionic currents remains hidden. In addition, we also assume for an experimental design which allows for the optional injection of external currents into the cell, thereby perturbing the action potential waveform. This opens the possibility for improvements in parameter identifiability, as the underlying model is required to simultaneously fit a larger set of independent transmembrane voltage data. Rather than directly fitting to the action potential data itself, we fit instead to the whole cell current required to generate each action potential data waveform. Our experience suggests that fitting directly to the action potential waveform is extremely difficult, particularly given the large numbers of parameters involved, due to the high number of local minima present. We find that the number of local minima may be reduced by modifying the original model equations so that the new model is effectively forced to follow the experimental data waveform. This is achieved through the use of a ‘data-clamp’ current injection, in which a simulation is performed in the presence of an additional compensatory current, icomp, added to the model equations for itot in (1), according to icomp ¼ gclamp ðVm
Vdata Þ

ð4Þ

where Vm is the model generated transmembrane voltage, Vdata is the experimentally observed transmembrane voltage at the same time instant, and gclamp is the fixed data-clamp conductance. The latter is typically set to a value of 300 mS/cm2, sufficiently high to ensure that the model waveform effectively follows (i.e. is ‘clamped’ to) the experimental data. Adding icomp into the model equations is equivalent to a proportional negative feedback system, ensuring that Vm effectively follows the Vdata waveform. The resulting itot in the presence of icomp then serves as the target whole cell current to be achieved by the model in the standard case (i.e. without icomp included in the equations). Least squares estimation of model parameters can then be carried out by minimising the sum of squares of the residual whole-cell current. Our experience is that the data-clamp approach is more effective in fitting ionic models to action potential data than by direct fits to membrane potential alone, due to the fact that icomp in (4) provides a direct measure of total membrane current discrepancy between model and data at the data action potential. This is because ionic currents typically depend on the transmembrane potential, as well as their own set of parameters. By forcing an ionic model to follow the data action potential waveform, a major source of variance between model and data has now been eliminated, since any difference in action potential waveforms will further alter the underlying ionic currents. It can be readily seen that the additional compensatory current required to be added to the data-clamped model must be a direct measure of total ionic current discrepancy between model and data given the data action potential waveform. Note that this membrane current discrepancy is not equivalent to simply comparing the difference between
Cm(dVdata/ dt) and itot in the unclamped case, as the action potential waveforms between model and data will, in general, be different. One advantage of this approach is that direct examination of this

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

412

Perturbed Transmembrane Potential 100

mV

50 0 -50 -100

Pseudo-Random Injection Current 40

inward µA/cm2

30 20 10 0 -10 -20

0

50

100

150

200

250

300

350

400

450

500

ms

Fig. 2. Effect of external pseudo-random current injection on transmembrane potential of BR model. Injected current is plotted according to the original BR convention for stimulus current: inward (depolarising stimuli) are shown as positive. Following the initial stimulus spike, a perturbing current is applied at mean hyperpolarising holding level of 5 mA/cm2. The top panel illustrates a series of complex ‘spike’ activation patterns in transmembrane potential, due to depolarising regions in the injected current waveform.

membrane current discrepancy may provide clues as to which ionic currents in the model need to have their parameters further modified. For example, examination of the reversal potential of this residual current, as well as its kinetics of activation and inactivation, may indicate which of the more dominant current components still need to be adjusted. If desired, the data-clamp protocol described above could also be carried out in the presence of external current injections, in order to sufficiently perturb the action potential waveform. When applied, the pseudo-random perturbation currents were delivered immediately following the excitatory stimulus. They were generated from a normal distribution with mean ilevel and standard deviation 2ilevel every 26.3 ms (i.e. 38 Hz), smoothed using cubic spline interpolation at all intermediate time values. Fig. 2 illustrates BR model behaviour when subjected to a pseudorandom current injection of ilevel ¼ 5 mA/cm2.

3. Parameter identifiability A model may be defined as globally identifiable in its parameters if given suitable experimental data, there exists a unique set of parameters for which the model is able to reproduce this data. In

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

413

the case of non-linear models, a priori determination of global identifiability is generally intractable, except in highly specific cases (Audoly et al., 2001). In contrast to the requirement of global identifiability, an alternative approach is that of local identifiability; defined as the existence of a unique set of parameters that satisfies the data within a local neighbourhood of parameter-space. According to least squares theory (Belsley, 1991), model parameters are locally identifiable provided the Hessian matrix is non-singular at the known solution. In the leastsquares case, the Hessian matrix, H, can be approximated from the Jacobian matrix, J, of the residual according to H ¼ 2JT J; provided the residual vector is sufficiently well approximated as linear with respect to the adjustable parameters. In this study, we examine the reciprocal condition value (1/l) of the Hessian. 1/l is defined as the ratio of the smallest to largest eigenvalues of H, and is equal to zero if the matrix is singular.

4. Least-squares optimisation by curvilinear gradient In non-linear least-squares optimisation, iterative techniques are employed which locally approximate the least-squares objective function q using a quadratic form qðxÞ ¼ c þ xT a þ 12xT Hx

ð5Þ

where x is an n  1 parameter vector, c is a scalar, a is the local n  1 gradient of the objective, and H is the local n  n Hessian matrix. In complex non-linear ODE systems models, including ionic models, a and H must be determined numerically at each step of the iteration. Direct numerical evaluation of H is generally not feasible, requiring O(n2) evaluations of the objective at each step. Instead H is usually approximated either from least-squares theory or by iterative updates using successive local gradient information, e.g. quasi-Newton methods (Fletcher, 1987). These latter methods only require an O(n) numerical evaluation of the gradient at each step. From leastsquares theory, a and H can be calculated using a ¼ 2JT r0 and H ¼ 2JT J; where J is the Jacobian matrix of the least squares residual with respect to the n parameters, and r0 is the residual vector evaluated using the current model parameters. Regardless of the method used to generate (5), two alternative courses of action can be taken to further advance to the minimum of the least squares objective. If (5) is assumed to be an exact representation of the objective, we can directly advance to the minimum using the full Newton step x ¼
H
1 a: Most optimisation methods employ a line search along this Newton direction in the hope that (5) eventually becomes sufficiently accurate to locate the minimum. However, if (5) is inexact far from the local x, then the Newton direction is of little value. Another difficulty arises when H is singular, as the Newton step is then undefined. An alternative is to search along the direction of local steepest-descent. This has the advantage that the objective will decrease locally at its maximum rate. Unfortunately, even if (5) is exact for all x, the method of steepest descent requires too many iterations in practice to descend to a minimum, since the steepest descent direction does not generally point toward this minimum. In contrast, only one iteration is required using the Newton direction when (5) is exact. Many investigators have suggested the possibility of combining the above two approaches by conducting a minimum search along a curvilinear path consisting initially of the steepest descent direction and terminating at the Newton step point. Such methods include the piece-wise linear

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

414

dogleg (Powell, 1970) and double dogleg (Dennis and Mei, 1979) trajectories as well as twodimensional subspace searches along parabolic arcs (Goldfarb, 1980). The popular Levenberg– Marquardt method for least-squares optimisation (More! , 1978) can be shown to be equivalent to such a curvilinear search which minimises (5) within a spherical trust region around the local x as the trust-region radius is increased (Fletcher, 1987). Another possibility is to conduct a curvilinear search along the steepest descent path of (5) all the way to the full-Newton step. Parameterising this trajectory as L(l): R1 -Rn ; where a is a scalar, it can be determined as the solution to the vector differential equation dLðaÞ ¼
rqðxÞ ¼
a
Hx da ¼
a
HLðaÞ

ð6Þ

the solution of which is x ¼ LðaÞ ¼ ðe
Ha
IÞH
1 a;

ð7Þ

where a is a scalar parameter ranging from 0 to N, defining a curvilinear trajectory L in Rn. Evaluation of (7) however, involves the calculation of a matrix exponential, which may be too expensive at each iteration step when compared with traditional methods, particularly for large H. Although such a gradient path approach has been previously suggested for the general minimisation problem (Botsaris and Jacobson, 1976; Vial and Zang, 1977), it has not been widely adopted, possibly due to the computational overhead involved. We have shown, however, that the method lends itself readily to the problem of least-squares minimisation, and that it surprisingly outperforms traditional least-squares approaches, particularly for large numbers of parameters (Dokos and Lovell, 2003). The following two relations can be readily shown to follow from Eq. (7): lim LðaÞ ¼
H
1 a

a-N

lim

a-0

ðfull Newton stepÞ;

dLðaÞ ¼
a ðsteepest-descent directionÞ: da

Thus, the trajectory given by (7) represents an optimal descent strategy combining the advantages of efficient local descent as well as off-local convergence to the quadratic minimum. Of all possible trajectories which interpolate between the steepest descent direction and the Gauss–Newton step, the steepest descent path described by Eq. (7) is optimal in the sense that it always follows the maximal descent direction, given that the quadratic approximation of Eq. (5) is accurate. Under far-from-local conditions, (5) may deviate from the real objective, however the optimal strategy is still to exploit the local quadratic approximation as far as possible, descending along its steepest descent trajectory until a minimum is found. In the Levenberg–Marquardt method, its curvilinear path also interpolates between the steepest descent direction and the Gauss-Newton step. However, its trajectory is equivalent to finding the minimum of the quadratic objective within a spherical trust region of given radius (Fletcher, 1987). This yields a different path which, in general, does not correspond to maximal descent direction at a given point, even if the underlying quadratic approximation is accurate.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

415

The steepest-descent trajectory described by Eq. (7) is illustrated in Fig. 3 for a quadratic form representing an arbitrary least-squares objective surface. Standard least-squares approaches compromise local descent efficiency by attempting to search for the objective minimum based on the far from local Newton step direction. Traditional steepestdecent methods, however, whilst maintaining local descent efficiency, compromise the far from local minimum direction as local conditions do not in-general reflect far-from-local behaviour. In order to conveniently calculate the curvilinear trajectory defined by (7), we can make use of the spectral decomposition of H, namely: H ¼ VDV
1 ; where D is a diagonal matrix whose elements, dii, are the eigenvalues of H, and V is the matrix whose columns are the corresponding eigenvectors. For a real symmetric matrix, D contains only real elements and V is orthogonal (i.e. VT ¼ V
1 ). It can be shown that (7) can be written in terms

Fig. 3. Curvilinear gradient path descent along a quadratic form representing an arbitrary least squares objective surface. The horizontal axes represent x1 and x2 adjustable parameters in simplified 2D parameter space. Also shown for comparison are the Newton and steepest-descent line directions. The curvilinear path represents the steepest descent trajectory along the quadratic surface. It terminates at the quadratic minimum (full Newton step), and is parallel to the steepest descent line-direction at the start of its descent.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

416

of this spectral decomposition as LðaÞ ¼ VMV
1 a
dii a

e

with M diagonal such that


1

; dii a0 dii ¼
a; dii ¼ 0

mii ¼

ð8Þ

Thus, the curvilinear trajectory given by (8) is analytically well-defined even when H is singular (dii ¼ 0). Note that since H is positive semi-definite in least-squares theory, then dii X 0. This means that (8) is always well-behaved, even for large a. At each step of the least-squares iteration, a minimum search along the curvilinear gradient defined by (8) was performed using Brent’s line-minimisation method (Press et al., 1994), modified to search along a curvilinear trajectory: the search being conducted on the extended positive reals from a ¼ 0 to N. A shortcoming of all descent methods of optimisation is that when a local minimum of the objective is reached, the method is not able to advance from this point, as there is no further downhill direction defined. To overcome this limitation, we implemented a global optimisation strategy through the use of a weighted least-squares residual scheme, and applied iterative reweighting when a local minimum was reached. In brief, if the local minimum yielded a non-zero objective function, the following weight update was applied: ¼ wðkÞ wðkþ1Þ i i þ sDwi ; ( "
"
"
# #  #) jri j
m1 2 jri j
m2 2 jri j
m3 2 þ exp þ exp ; Dwi ¼ exp s s s

ð9Þ

ðkþ1Þ where wðkÞ are the weights of the ith residual, ri ; before and after the kth update i ; wi respectively, and m1, m2, m3, s and s are positive scalars defined according to

m1 ¼ r1 minðjrjmax ; 2jrjmedian Þ; m2 ¼ r2 minðjrjmax ; 2jrjmean Þ; m3 ¼ ð0:2r3 þ 0:8Þjrjmax ; s ¼ jrjSD ; "rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # P ðkÞ 2 5 min 10 wi ri ; 10 s¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ðDwi ri Þ2 þ 10
12

;

ð10Þ

where |r|max, |r|mean and |r|median denote respectively the maximum, mean, and median values of the absolute residual, r1, r2, r3 are random variables from a uniform distribution in [0,1], and s is calculated such that the new objective evaluated with the updated weight is 10-fold larger than the previous objective and has a ceiling value of 1010 (to prevent scaling problems associated with extremely large weights). The re-weighting scheme described in Eqs. (9) and (10) was chosen so that absolute residuals corresponding, on average, to the median, the mean, and 90% of the maximum of all the absolute residuals are weighted the most heavily. This is achieved by summing

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

417

three bell-shaped weight functions, with some degree of randomness in their generation so that repeating cycles of weights are not endlessly applied. Maximally weighting the mean, median and 90% maximum absolute residuals addressed the case of fitted residuals which may have abnormally large maximal values compared to their mean or median. In such a case, a large weight update would be applied to possibly only very few data points at the expense of remaining portions of the residual. Emphasising mean and median values to the weight function ensures that a wider set of residual points is weighted. After a few iterations, k, of the updating scheme, the weight is reset to its initial value. This reset occurs whenever the new weighted objective value decreases by more than 70% of its value at the start of the current iteration, or a maximum of 10 updates have already been applied. If after resetting the weight, and following another search, the current minimum objective value has not been improved on, the parameters are randomly perturbed around the best objective minimum parameter set found, and a new search is initiated. The initial weight is given by 10 : ð11Þ wð0Þ i ¼ 10 þ jri j This weighting function applies an initial weight of unity to small values of the absolute residual, and a smaller weight to larger values of the initial absolute residual, since these values may contribute an inordinate amount to the overall objective. Finally, all parameters were constrained to lie within physiological ranges by restricting their value to lie between 7100% of their default value given in the appendix. These constraints were implemented using a transformation approach (Fletcher, 1987), in which all 63 parameters, pi ; were transformed into new unconstrained parameters, xi, using the transformation equations

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi
Li
1 ; xi ¼ sin Ui
L i pi ¼ Li þ ðUi
Li Þ sin2 ðxi Þ;

ð12Þ

where Ui ; Li denote the upper/lower limit for model parameter i (set, respectively, at 7100% of default). The curvilinear gradient optimisation method described above was performed on the transformed parameter set xi. All computations were carried out on a Pentium PC using Matlab software (v6.5, The Mathworks, inc).

5. Results 5.1. Effect of parameter variations In order to assess the effect of modest parameter perturbations on action potential wave-shape, all 63 parameters of the BR model (i.e. pi ) were randomly perturbed from their default values using a uniform probability distribution centred around the default value of each parameter. Two sets of simulations were run, consisting of 10 runs of randomly perturbed parameters. In the first set of simulations, the maximum range of the parameter deviations was set to 710% of each default value, corresponding to a standard deviation of 5.8%. In the second set of simulations, parameters were perturbed by 730% of default, corresponding to a standard deviation of 17.3%.

ARTICLE IN PRESS 418

S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431 40 10% max 20

E (mV)

0 -20 -40 -60 -80 -100

0

50

100

150

200

250

300

350

400

450

500

400

450

500

time (ms) 40 30% max

20

E (mV)

0 -20 -40 -60 -80 -100

0

50

100

150

200

250

300

350

time (ms)

Fig. 4. Effect of modest parameter variations on BR action potential. For both panels, the standard BR action potential generated from default parameter values is shown in bold. In each panel, ten sets of 63 parameters were generated from a uniform random distribution centred around the default values. Maximum parameter deviations were set at 710% (upper panel) and 730% (lower panel) from default.

Results for both runs are shown in Fig. 4, where it is evident that at 710% maximum deviation, the action potential has somewhat preserved its characteristic shape. At maximum deviations of 730% however, characteristic ventricular action potential wave-shape was in general not maintained. This suggests that modest parameter uncertainties of up to 730% in cardiac ionic models may lead to substantial inaccuracies in model behaviour, particularly when large numbers of parameters are present. 5.2. Parameter identification Table 1 illustrates the value of the reciprocal condition value (1/l) of the Hessian at the known solution point for all the parameters of the BR model in the case of n ¼ 1; 2, and 3 total

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

419

Table 1 Reciprocal condition (1=l) value for four experimental designs described in text Experimental design

n

1=l

Default action potential (AP) Default AP+1 perturbed waveform Default AP+2 perturbed waveforms

1 2 3

3.6310e-14 5.1761e-11 1.3691e-10

To generate these values, the data-clamp objective function was used at the default BR parameter solution point. Small 1=l values imply a near-singular Hessian matrix, therefore locally unidentifiable parameters at the known solution.

experimental waveforms simultaneously fitted. These correspond, respectively, to the cases of a single action-potential waveform (i.e. no perturbing current), one and two additional pseudorandom perturbing current injections. For these sets of waveforms, the mean level of perturbing current, ilevel, was adjusted for the ith perturbation run (i ¼ 1; y; n; where n is the total number of waveforms), according to i
1 ; mA=cm2 ilevel ¼ 10 n
1 which corresponded to mean levels of hyperpolarising injections. It can be readily seen that 1/l remains sufficiently small for each additional current injection experiment, that the Hessian is effectively singular in each case. This suggests that the BR model is over-determined in all 63 parameters, and that additional data from perturbed action potentials is unlikely to improve this identifiability problem. Fig. 5 illustrates the result of one parameter fit to a single action potential waveform generated using the default BR model. In this fit (termed ‘‘Fit 1’’), all 63 model parameters were randomly assigned using a uniform probability distribution centred around the default value of each parameter, with a total range of 730% of default. As can be seen from Fig. 5, the action potential was well fit by the optimisation techniques presented earlier, exhibiting a fitted root-mean square (rms) error of only 0.37 mV. All 63 parameter values obtained from Fit 1 are shown in Table 2, with a summary of parameter deviations from default given in Table 3. The results of this fit reveal that despite obtaining a good fit to the original action potential waveform, fitted parameters were substantially different from known values of the default model (mean deviation, 14.5%, see Table 3). In addition, all ionic currents as well as intracellular calcium concentration were not able to be properly reconstructed, as shown in the lower panels of Fig. 5, even though the underlying action potential was well fit. Fig. 6 illustrates the results of an additional fit to the BR model, termed ‘‘Fit 2’’. In this case, the initial model parameter values were selected randomly from a uniform distribution ranging between 710% of default values. Although these initial parameters represent modest excursions from the default BR model, the resulting fit to the single action potential data (Fig. 6, top panel) did not exhibit accurate reconstructions of underlying ionic mechanisms (Fig. 6, lower panels), with a mean parameter deviation of 8.4% from known values (Table 3). Results from Fits 1 and 2 indicate that good fits to single action potential data alone is a poor indicator of the accuracy of underlying ionic current reconstructions, particularly given the over-determined nature of the model.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

420

40 20

default model fit 1

E

0 RMS error: 0.36852 mV mV

-20 -40 -60 -80 -100 100 istim

µA/cm

2

0

-100

iNa

-200 10

µA/cm

2

5

iK1 ix1

0

is -5 8 [Ca]i

µM

6 4 2 0 0

50

100

150

200

250

300

350

400

450

500

ms

Fig. 5. ‘‘Fit 1’’ results of BR model least-squares fit to single action potential data (top panel). All 63 BR model parameters were initially assigned random values from a uniform distribution ranging between 730% of parameter defaults. The fit reveals that default BR action potential waveforms can be generated even though underlying membrane currents and ionic concentrations are not properly reconstructed. In all panels, fitted model behaviour (dashed lines) is compared to the default BR model (solid lines). A list of all parameter values obtained for this fit is shown in Table 2, with summary of comparison to known model parameters shown in Table 3.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

421

Table 2 List of all 63 fitted parameter values obtained for three fits to default BR action potential data BR Parameter

Default Value

Fit 1

Fit 2

Fit 3

P1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p14 p15 p16 p17 p18 p19 p20 p21 p22 p23 p24 p25 p26 p27 p28 p29 p30 p31 p32 p33 p34 p35 p36 p37 p38 p39 p40 p41 p42 p43 p44 p45 p46 p47

1.4000 0.0400 85.0000 0.0800 53.0000 0.0400 0.0700 23.0000
0.0400 0.1973 0.0400 77.0000 0.0400 4.0000 0.0030 50.0000 0.0900
82.3000
13.0287
1.0000e-7 0.0700 1.0000e-7 0.0018 0.0830 0.0570 0.0578 8.7142e-4
0.0600
0.0400 2.2255
1.0000 47.0000
0.1000 0.7096
0.0560 5.4980e-10
0.2500 10.7578
0.0820 6.3281 0.0011
0.2500
0.2000 5.9565e+6 7.3598
0.1000 24.5325

1.2682 0.0301 80.1449 0.0862 42.1964 0.0344 0.0759 24.7872
0.0300 0.1373 0.0406 53.4430 0.0289 3.4270 0.0021 48.0127 0.0804
101.2821
15.2003
8.6332e-8 0.0762 9.8475e-8 0.0018 0.0561 0.0557 0.0362 9.9572e-4
0.0682
0.0514 2.3028
1.1819 35.2745
0.0768 0.6222
0.0549 3.7419e-10
0.2277 11.7484
0.1159 7.9550 9.1110e-4
0.2503
0.1660 5.7503e+6 6.8729
0.1002 20.1666

1.3241 0.0351 84.2431 0.0927 76.9697 0.0541 0.0697 22.9082
0.0403 0.2173 0.0395 82.0338 0.0416 4.2335 0.0027 49.6158 0.0840 -70.4970 -12.1706 -1.0113e-7 0.0735 8.5349e-8 0.0020 0.0648 0.0692 0.0485 8.8765e-4
0.0494
0.0375 2.2737
1.0207 46.7906
0.0989 0.7727
0.0555 5.2369e-10
0.2416 11.1054
0.0921 6.6709 0.0012
0.2282
0.2332 5.9800e+6 6.9464
0.0952 24.7577

1.6193 0.0496 85.5781 0.0809 61.8025 0.0427 0.0656 26.7826
0.0485 0.2077 0.0390 77.3874 0.0379 3.6784 0.0028 51.9292 0.0964
101.8665
14.0114
9.3329e-8 0.0665 9.1069e-8 0.0021 0.0759 0.0533 0.0679 6.8649e-4
0.0688
0.0492 2.4388
1.3564 41.5389
0.1052 0.7256
0.0535 5.5483e-10
0.2443 10.9365
0.0821 6.2252 8.8093e-4
0.2660
0.2165 5.7362e+6 6.3447
0.1048 21.8345

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

422

Table 2 (continued) BR Parameter

Default Value

p48 p49 p50 p51 p52 p53 p54 p55 p56 p57 p58 p59 p60 p61 p62 p63

0.0697
0.0100
0.0720 0.6977 0.0037
0.0170 0.0500 0.1108 1.4383e-4
0.0080 0.1500 0.0150 1.4391
0.0200
0.2000 403.4300

Fit 1 0.0672
0.0065
0.0766 0.6189 0.0046
0.0154 0.0454 0.0965 1.4770e-4
0.0086 0.1530 0.0179 1.5505
0.0217
0.2290 351.4584

Fit 2 0.0640
0.0087
0.0909 0.6978 0.0041
0.0174 0.0550 0.1115 1.6182e-4
0.0082 0.1812 0.0175 1.5541
0.0196
0.2002 440.7546

Fit 3 0.0758
0.0124
0.0862 0.7770 0.0041
0.0159 0.0573 0.1263 1.5021e-4
0.0094 0.1796 0.0186 1.3451
0.0213
0.1741 424.7257

Default values for all parameters, equivalent those values of the standard BR formulation, are shown in column 2.

Table 3 Summary of fitted parameter data for all three fits to the BR model described in the text

Fit 1 Fit 2 Fit 3

Maximum absolute deviation in all 63 fitted parameters (relative to default values) (%)

Mean absolute deviation in all 63 fitted parameters (relative to default values) (%)

41.3 45.2 35.6

14.5 8.4 10.5


 pi
p  Maximum absolute deviation relative to known parameter defaults was calculated using maxi   i   100%; pi where p is the ith fitted parameter value and p the known (i.e. default) value for that parameter. Note that all i

i

maximum deviations from default were well within the 100% upper/lower range constraints of parameters, implemented using the transform equation (12). Mean absolute relative deviation was calculated in a similar fashion.

5.3. Reconstruction of hidden ionic currents In order to establish whether underlying ionic mechanisms have, in fact, been accurately reconstructed when fitting to transmembrane potential data, the experimenter has the option of perturbing the whole-cell recording by means of external current injections. The additional perturbations may act as predictive tests, verifying whether these parameters are able to reproduce the additional data. If the parameters exhibit an unsatisfactory fit to the new data, then this data may be included in subsequent fitting processes, in order to refine the accuracy of the underlying ionic current reconstructions. An example of such a fit, termed ‘‘Fit 3’’, is shown in Fig. 7. Similar

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

423

40 default model fit 2

E

20 0

RMS error: 0.30887 mV mV

-20 -40 -60 -80 -100 100 istim

µA/cm

2

0

-100

iNa

-200 5

µA/cm

2

iK1 ix1

0

is -5 8 [Ca] i

µM

6 4 2 0

0

50

100

150

200

250 ms

300

350

400

450

500

Fig. 6. ‘‘Fit 2’’ results of BR model least-squares fit to single action potential data (top panel). All 63 BR model parameters were initially assigned random values from a uniform distribution ranging between 710% of parameter defaults. The fit reveals that default BR action potential waveforms can be generated even though underlying membrane currents and ionic concentrations are not properly reconstructed. In all panels, fitted model behaviour (dashed lines) is compared to the default BR model (solid lines). A list of all parameter values obtained for this fit is shown in Table 2, with summary of comparison to known model parameters shown in Table 3.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

424

40 20

default model fit 3

E

0 RMS error: 0.79275 mV mV

-20 -40 -60 -80 -100 100 istim

µA/cm

2

0

-100

iNa

-200 5

µA/cm

2

iK1

ix1

0

is -5 8 [Ca] i

µM

6 4 2 0

0

50

100

150

200

250 ms

300

350

400

450

500

Fig. 7. ‘‘Fit 3’’ results of BR model least-squares fit to two sets of simultaneous data: a single action potential and perturbed waveform generated using a mean hyperpolarising injected current of 10 mA/cm2. The fit reveals that the default BR action potential waveform can be reproduced by this model (top panel), as well as adequate reconstructions of default BR membrane currents and ionic concentrations (lower panels). All 63 BR model parameters were initially assigned random values from a uniform distribution ranging between 710% of parameter defaults. In all panels, fitted model behaviour (dashed lines) is compared to the default BR model (solid lines). A list of all parameter values obtained for this fit is shown in Table 2, with summary of comparison to known model parameters shown in Table 3.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

425

to Fit 2 of Fig. 6, initial model parameters where randomly generated from a uniform distribution within 710% of default values. In this case, however, simultaneous fits were obtained to two sets of data, a default action potential as well as a perturbed transmembrane potential waveform generated with a mean hyperpolarising current level of 10 mA/cm2. This time, the resulting fitted Perturbed Transmembrane Potential default model fit 1 fit 2 fit 3

200

mV

100 0 -100 0

50

100

150

200

250 ms

300

350

400

450

500

10

µA/cm

2

0 Fully-activated ix1

-10 -20 -30 -120 6

x 10

-100

-80

-60

-3

-40 mV

α

x1

0

20

40

-20

0

20

40

-20

0

20

40

rate

α x1

4

-20

2

0 -120

-100

-80

-60

-40 mV

0.01

β x1

0.008

β x1 rate

0.006 0.004 0.002 0 -120

-100

-80

-60

-40 mV

Fig. 8. Accuracy of ionic current reconstructions for all three fits obtained. Top panel illustrates a perturbation simulation in which a pseudo-random current with a mean hyperpolarizing level of 10 mA/cm2 is injected into the model. The lower three panels illustrates the reconstructed fully activated current-voltage relation for ix1, and the reconstructed ax1 ; bx1 rate coefficients of the x1 gating variable for all three fits (dashed lines) as well as the default model (solid lines).

ARTICLE IN PRESS 426

S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

model depicted accurate reconstructions of hidden ionic mechanisms, as shown in the lower panels of Fig. 7. Despite the improved reconstruction in ion currents, however, fitted parameters still exhibited a mean deviation of 10.5% from known values (Table 3). This suggests that this experimental design is still not able to uniquely identify model parameters, as evidenced by the near-singular Hessian matrix (1=l ¼ 5:2e
11) for the two-waveform experiment (Table 1). To examine whether current-injection perturbations are able to provide an indication of the accuracy of underlying ionic reconstructions, Fig. 8 (top panel) illustrates the results of a pseudorandom current perturbation applied to the BR model using all three fits from Figs. 5–7, as well as the default BR model itself. It can be seen from this result that in contrast to Fit 3, Fits 1 and 2 parameters were not able to reproduce the perturbation data, since these fits yielded inaccurate reconstruction of hidden ionic mechanisms (see Figs. 5 and 6). In order to assess the accuracy of the underlying current kinetics and fully activated amplitudes obtained from all three fits, the lower panels of Fig. 8 illustrate data obtained for the timedependent ix1 current. It can be seen that only Fit 3 was able to accurately reconstruct the fully activated ix1 current–voltage characteristic, as well as the associated ax1 ; bx1 voltage characteristics of the rate coefficients of the x1 gating variable. These characteristics were closely reproduced, despite the fact that underlying parameters were different from known default values. In the case of Fit 3, there was a deviation in the value of the bx1 rate coefficient from default at large negative membrane voltages below –60 mV (lower panel, Fig. 8). This inaccuracy was also evident on close examination of the perturbed action potential data (top panel), where large negative voltages in the transmembrane potential could not accurately be fitted by these sets of parameters, due to the presence of multiple local minima in the least-squares objective. In such cases, it may be possible to design experiments which inject perturbation currents to elicit transmembrane responses in the voltage region known to be inaccurately fitted.

6. Discussion The results of this study suggest that the BR model is an over-determined system in all 63 parameters defined, when fitting to a single action potential waveform. It also suggests that moderate parameter variations of up to 730%, well within the range of typical biological variability, may lead to inaccurate reconstructions of the cardiac action potential (see Fig. 4, lower panel). It is therefore highly desirable to implement methods which can accurately identify model parameters against data obtained from specific cells. This is, however, confounded by the overdetermined nature of the BR model. Local parameter identifiability in the BR model could, however, be moderately improved by the addition of pseudo-random current injection data. Despite the additional data, the model remained essentially indeterminate. Since the BR model represents a simple example of a cardiac ionic model, especially when compared to more recent models which employ a host of ion channels, pumps and exchangers (Luo and Rudy, 1994), metabolic processes (Ch’en et al., 1998), ATP energetics (Michailova and McCulloch, 2001), as well as excitation–contraction coupling (Jafri et al., 1998), it may be reasonably asserted that such models will be over-determined in their full set of parameters. It may, however, still be possible to accurately reconstruct

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

427

underlying ionic current kinetics and amplitudes through the use of a simple experimental design consisting of external current injections to perturb the action potential waveform. In this study, we have shown (see Fig. 8) that such perturbation data is able to qualitatively reveal deficiencies in BR model parameters that have been previously fitted to single action potential data alone. Thus, although it is not possible to uniquely identify all parameters in an ionic model, what is important to the experimentalist is the reconstruction of hidden ionic current amplitudes and kinetics which underlie the action potential waveform. We report that the ability of a model to reproduce additional perturbed current injection waveforms is a good indicator of how well underlying currents have been reconstructed. This is despite the fact that individual model parameters have not been uniquely identified, due to its over-determined nature. Through the use of a custom least-squares optimisation routine described in this study, fits to single action potential data could be readily obtained (Fits 1 and 2). Table 2 illustrates three additional distinct sets of 63 model parameters which were able to accurately reproduce the default BR action potential waveform. In terms of fitting model parameters to multiple sets of action potential data, however, we find in practice that simultaneous fitting presents significant difficulties in terms of computational efficiency. That is, global optimisation of model parameters is much more computationally onerous using multiple data, due to the fact that the number of local minima in parameter-space increases substantially with each additional data set. One approach, which may be feasible in practice, is to simply fit two simultaneous sets of data: the default action potential data and one additional perturbed waveform, as was the case with Fit 3. Results of this study suggest that such an approach is capable of yielding accurate reconstructions of ionic mechanisms underlying the cardiac action potential.

7. Conclusions We have presented a novel experimental design that allows for the computational reconstruction of underlying ionic mechanisms in essentially over-determined cardiac ionic models. The kinetics and magnitudes of hidden ionic currents can be recovered through the use of suitable pseudo-random current injections. It is anticipated that this method may be applied to experimental data from real myocytes to generate cell-specific parameter fits to a host of cardiac cell types. Appendix The equations defining the complete BR ionic model of the ventricular action potential, together with the specific parameterisation adopted in this study, are given below. A total of 63 adjustable parameters were employed, denoted by pi (i = 1, 2, y63), with default values given. In some cases, the equations were reformulated (by redefining coefficients) in order to remove any obvious redundancies in parameters (see text). All equations shown below are mathematically equivalent to the original BR descriptions.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

428

A.1. Ion currents iK1 ¼

p1 ½exp½p2 ðVm þ p3 Þ
1 exp½p4 ðVm þ p5 Þ þ exp½p6 ðVm þ p5 Þ p7 ðVm þ p8 Þ þ ; 1
exp½p9 ðVm þ p8 Þ

ðA:1Þ

where p1 ¼ 1:4 mA/cm2, p2 ¼ 0:04 mV
1, p3 ¼ 85 mV, p4 ¼ 0:08 mV
1, p5 ¼ 53 mV, p6 ¼ 0:04 mV
1, p7 ¼ 0:07 mA/cm2, p8 ¼ 23 mV, and p9 ¼ 0:04 mV
1. ix1 ¼

x1 p10 ½exp½p11 ðVm þ p12 Þ
1 ; exp½p13 Vm 

ðA:2Þ

where p10 ¼ 0:1973 mA/cm2, p11 ¼ 0:04 mV
1, p12 ¼ 77 mV, and p13 ¼ 0:04 mV
1. iNa ¼ ðp14 m3 hj þ p15 ÞðVm
p16 Þ;

ðA:3Þ

where p14 ¼ 4 mS/cm2, p15 ¼ 0:003 mS/cm2, and p16 ¼ 50 mV. is ¼ p17 df ðVm
p18
p19 ln½Cai Þ;

ðA:4Þ

where p17 ¼ 0:09 mS/cm2, p18 ¼
82:3 mV, and p19 ¼
13:0287 mV/ln(M). A.2. Membrane potential
 dVm 1 ðiK1 þ ix1 þ iNa þ is
istim Þ; ¼
Cm dt

ðA:5Þ

where Cm is held fixed at 1 mF/cm2. A.3. Intracellular calcium concentration d½Cai ¼ p20 is þ p21 ðp22
½Cai Þ; dt where p20 ¼
1e
7 M cm2/nC, p21 ¼ 0:07 ms
1, and p22 ¼ 1e
7 M.

ðA:6Þ

A.4. Gating variables dx1 ¼ ax1 ð1
x1 Þ
bx1 x1 ; dt

ðA:7Þ

dm ¼ am ð1
mÞ
bm m; dt

ðA:8Þ

dh ¼ ah ð1
hÞ
bh h; dt

ðA:9Þ

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

429

dj ¼ aj ð1
jÞ
bj j; dt

ðA:10Þ

d ½d  ¼ ad ð1
dÞ
bd d; dt

ðA:11Þ

df ¼ af ð1
f Þ
bf f ; dt

ðA:12Þ

with all a, b rate coefficients (in units of ms
1) given as: p23 exp½p24 Vm  ax1 ¼ ; exp½p25 Vm  þ p26

ðA:13Þ

(p23 ¼ 0:0018 ms
1, p24 ¼ 0:083 mV
1, p25 ¼ 0:057 mV
1, and p26 ¼ 0:0578) bx1 ¼

p27 exp½p28 Vm  ; exp½p29 Vm  þ p30

ðA:14Þ

(p27 = 8.7142e
4 ms
1, p28 =
0.06 mV
1, p29=
0.04 mV
1, and p30=2.2255) am ¼

p31 ðVm þ p32 Þ ; exp½p33 ðVm þ p32 Þ
1

ðA:15Þ

(p31 ¼
1 mV
1 ms
1, p32 ¼ 47 mV, p33 ¼
0:1 mV
1) bm ¼ p34 exp½p35 Vm ;

ðA:16Þ

(p34 ¼ 0:7096 ms
1, p35 ¼
0:056 mV
1) ah ¼ p36 exp½p37 Vm ; (p36 ¼ 5:498e
10 ms
1, p37 ¼
0:25 mV
1) p38 bh ¼ exp½p39 Vm  þ p40

ðA:17Þ

ðA:18Þ

(p38 ¼ 10:7578 ms
1, p39 ¼
0:082 mV
1, and p40 ¼ 6:3281) aj ¼

p41 exp½p42 Vm  ; exp½p43 Vm  þ p44

(p41 ¼ 0:0011 ms
1, p42 ¼
0:25 mV
1, p43 ¼
0:2 mV
1, and p44 ¼ 5:9565e6) p45 bj ¼ exp½p46 Vm  þ p47

ðA:19Þ

ðA:20Þ

(p45 ¼ 7:3598 ms
1, p46 ¼
0:1 mV
1, and p47 ¼ 24:5325) ad ¼

p48 exp½p49 Vm  ; exp½p50 Vm  þ p51

ðA:21Þ

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

430

(p48 ¼ 0:0697 ms
1, p49 ¼
0:01 mV
1, p50 ¼
0:072 mV
1, and p51 ¼ 0:6977) bd ¼

p52 exp½p53 Vm  ; exp½p54 Vm  þ p55

ðA:22Þ

(p52 ¼ 0:0037 ms
1, p53 ¼
0:017 mV
1, p54 ¼ 0:05 mV
1, and p55 ¼ 0:1108) af ¼

p56 exp½p57 Vm  ; exp½p58 Vm  þ p59

ðA:23Þ

(p56 ¼ 1:4383e
4 ms
1, p57 ¼
0:008 mV
1, p58 ¼ 0:15 mV
1, and p59 ¼ 0:015) bf ¼

p60 exp½p61 Vm  ; exp½p62 Vm  þ p63

ðA:24Þ

(p60 ¼ 1:4391 ms
1, p61 ¼
0:02 mV
1, p62 ¼
0:2 mV
1, and p63 ¼ 403:43) A.5. Initial conditions Vm ¼
83:3 mV, [Ca]i=1.86e
7 M, x1 ¼ 0:1644; m ¼ 0:01; h ¼ 0:9814; j ¼ 0:9673; d ¼ 0:0033; f ¼ 0:9884:

References Audoly, S., Bellu, G., D’Anglio, L., Saccomani, M.P., Cobelli, C., 2001. Global identifiability of nonlinear models of biological systems. IEEE Trans. Biomed. Eng. 48 (1), 55–65. Beeler, G.W., Reuter, H., 1977. Reconstruction of the action potential of ventricular myocardial fibres. J. Physiol. 268, 177–210. Belsley, D.A., 1991. Conditioning Diagnostics: Collinearity and Weak Data in Regression. Wiley, New York. Botsaris, C.A., Jacobson, D.H., 1976. A newton-type curvilinear search method for optimization. J. Math. Anal. Appl. 54, 217–229. Ch’en, F.F.T., Vaughan-Jones, R.D., Clarke, K., Noble, D., 1998. Modelling myocardial ischemia and reperfusion. Prog. Biophys. Mol. Biol. 69, 515–538. Dennis, J.E., Mei, H.H.W., 1979. Two new unconstrained optimization algorithms which use function and gradient values. J. Optim. Theory Appl. 28, 453–482. Dokos, S., Lovell, N.H., 2003. A curvilinear gradient path method for optimization of biological systems models. In: Feng, D.D., Carson, E.R. (Eds.), Proceedings of the fifth IFAC Symposium on Modelling and Control in Biomedical Systems 2003. Elsevier, Oxford, pp. 203–208. Feng, J., Yue, L., Wang, Z., Nattel, S., 1998. Ionic mechanisms of regional action potential heterogeneity in the canine right atrium. Circ. Res. 83, 541–551. Fletcher, R., 1987. Practical Methods of Optimization, 2nd Edition.. John Wiley, Chichester. Goldfarb, D., 1980. Curvilinear path steplength algorithms for minimization which use directions of negative curvature. Math. Programming 18, 31–40. Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Phys. 117, 500–544. Hunter, P.J., Kohl, P., Noble, D., 2001. Integrative models of the heart: achievements and limitations. Philos. Trans. Roy. Soc. London, Ser. A 359, 1049–1054. Jafri, M.S., Rice, J.J., Winslow, R.L., 1998. Cardiac Ca2+ dynamics: The roles of ryanodine receptor adaptation and sarcoplasmic reticulum load. Biophys. J. 74, 1149–1168.

ARTICLE IN PRESS S. Dokos, N.H. Lovell / Progress in Biophysics & Molecular Biology 85 (2004) 407–431

431

Liu, D., Gintant, G.A., Antzelevitch, C., 1993. Ionic basis for electrophysiological distinctions among epicardial, midmyocardial, and endocardial myocytes from the free wall of the canine left ventricle. Circ. Res. 72, 671–687. Luo, C.H., Rudy, Y., 1994. A dynamical model of the cardiac ventricular action potential. I. simulations of ionic currents and concentration changes. Circ. Res. 74, 1071–1096. Michailova, A., McCulloch, A., 2001. Model study of ATP and ADP buffering, transport of Ca2+ and Mg2+, and regulation of ion pumps in ventricular myocyte. Biophys. J. 81, 614–629. Mor!e, J.J., 1978. The Levenberg–Marquardt algorithm: implementation and theory. In: Watson, G.A. (Ed.), Lecture Notes in Mathematics, Vol. 630. Springer, Berlin, pp. 105–116. Noble, D., Bett, G.C.L., 1993. Reconstructing the heart: a challenge for integrative physiology. Cardiovas. Res. 27, 1701–1712. Powell, M.J.D., 1970. A hybrid method for nonlinear equations. In: Rabinowitz, P. (Ed.), Numerical Methods for Nonlinear Algebraic Equations. Gordon and Breach, London, pp. 87–114. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1994. Numerical Recipes in C: the art of scientific computing, 2nd Edition.. Cambridge University Press, New York. Vial, J.-P., Zang, I., 1977. Unconstrained optimisation by approximation of the gradient path. Math. Oper. Res. 2, 253–265.