Abstract

Markov models (MMs) represent a generalization of Hodgkin–Huxley models. They provide a versatile structure for modelling single channel data, gating currents, state-dependent drug interaction data, exchanger and pump dynamics, etc. This paper uses examples from cardiac electrophysiology to discuss aspects related to parameter estimation. (i) Parameter unidentifiability (found in 9 out of 13 of the considered models) results in an inability to determine the correct layout of a model, contradicting the idea that model structure and parameters provide insights into underlying molecular processes. (ii) The information content of experimental voltage step clamp data is discussed, and a short but sufficient protocol for parameter estimation is presented. (iii) MMs have been associated with high computational cost (owing to their large number of state variables), presenting an obstacle for multicellular whole organ simulations as well as parameter estimation. It is shown that the stiffness of models increases computation time more than the number of states. (iv) Algorithms and software programs are provided for steady-state analysis, analytical solutions for voltage steps and numerical derivation of parameter identifiability. The results provide a new standard for ion channel modelling to further the automation of model development, the validation process and the predictive power of these models.

1. Introduction

Almost immediately following Hodgkin and Huxley's (HH) paper (Hodgkin & Huxley 1952), Markov models (MMs) were used to model ligand-binding ion channels (e.g. Armstrong 1971; Colquhoun 1973), but only recently have they been applied in whole-cell models (e.g. Winslow et al. 1999). As MMs are a generalization of HH models, they are more versatile and can represent state transitions that are dependent on the current state of the channel, thereby being able to match the experimental observations on gating currents (Varghese & Boland 2002). State-dependent drug interactions can obviously be implemented using MMs (e.g. Armstrong 1971), but are also possible using HH formulations (Starmer et al. 1984). The latest review on MMs (Rudy & Silva 2006) has emphasized these advantages of MMs, but has not covered the numerical issues with parameter estimation and model layout, which we discuss in this paper.

This paper focuses on voltage-gated ion channels for which there are only very limited biophysical data available compared with transporters and exchangers. For the latter, information exists for the estimation of model structure and a number of parameters, such as binding and dissociation constants of ions (e.g. Terkildsen et al. 2007).

Fredkin et al. (1985) have shown that ideal single channel data obtained from voltage-clamp steps provide all information necessary for fitting voltage-gated ion channel models, but with real, noisy data the analysis of first latencies is extremely difficult (Magleby 1992). As whole-cell current measurements are widely used to describe ion channel kinetics under physiological conditions, we concentrate on these for our investigations.

An algorithm for parameter estimation based on the log-likelihood method has been proposed by Milescu et al. (2005), who made their software publicly available. However, our own experience (Fink et al. 2008b) and remarks in other publications on using a ‘manual tweak’ (Vecchietti et al. 2006) indicate that parameter estimation is still a crucial problem when developing MMs (at least for voltage-gated ion channels). Problems with parameter estimation (independent of the algorithm used) make it impossible to determine the correct structure or layout of a MM (and there can also exist a priori indistinguishable models; see Kienker 1989). Therefore, the connection to the protein structure is not an emerging property of the modelling, but dependent on the modeller's input.

Last but not least, MMs have been associated with high computational cost owing to their large number of states (and therefore large number of differential equations), which presents an obstacle for both parameter estimation and multicellular whole-organ simulations (Hunter et al. 2003). The question arises about how we can make best use of the current knowledge and experimental data to support multiscale modelling.

In this paper, we examine various Markov and HH models of voltage-gated ion channels regarding the identifiability of their parameters, i.e. the possibility of finding a unique set of parameters which provides the best fit of the model to experimental data. The dependency of the identifiability on the available experimental data is discussed, and we present a shortened voltage-step clamp protocol to improve the parameter estimation process and provide the possibility of model validation. Furthermore, we discuss advantages and disadvantages of using HH formulations and of reducing the number of states of MMs.

2. Hodgkin–Huxley and Markov model formalisms

MMs and HH models represent a way to describe ion channel kinetics. The basis for the models is that the ion channel proteins have different confirmational states (open, closed and inactivated) that can be found when investigating the molecular structure (Isacoff et al. 1993; Doyle 2004). The models can have multiple open, closed and/or inactivated states, but only the open states are conducting states. The structure of the MM can be described by the states and the transitions between these states. The models in figure 1 denote various MMs, indicating the versatility of the combination of states and transitions. However, the model in figure 1a is a special case of a MM, having a rectangular grid of states with very particular transition rates assigned. This is the representation of an HH model using the MM formalism.

The structure of four different model types using the MM formalism. (a) The layout of HH-type models such as INa by Hodgkin & Huxley (1952) or Ito by Campbell et al. (1993b)—note the small number of rate constants. (b–d) MM formulations of INa (Clancy & Rudy 2002), Ito (Campbell et al. 1993b) and ICaL (Faber et al. 2007). MMs can have irregular shapes and the transition rates can be selected independently (observing microscopic reversibility). We re-use the structure of the model in (b) for our own (identifiable) INa model presented in this paper. C, closed; I, inactivated; CI, closed inactivated; IV, voltage inactivated; O, open states. The α, β, γ, δ, λ, Φ, ϕ and ω denote the respective transition rates between states.

The main advantage of the MM versus the HH formalism is the large degree of freedom in the model structure that enables it to replicate measurements of gating currents (Armstrong & Bezanilla 1977; Varghese & Boland 2002), which show a state dependence of activation and inactivation, contradicting the assumption of independence of gating variables in the HH approach. However, the importance of this effect on whole-cell currents has not yet been demonstrated.

3. On the determination of the structure of Markov-type models

The higher flexibility in the layout of MMs comes with the issue of finding the optimal model structure/layout for each ion channel. There have been a number of publications proposing optimal layouts of MMs for specific ion channels and their advantages over HH models (e.g. Horn & Vandenberg 1984; Chay 1991; Vandenberg & Bezanilla 1991; Destexhe et al. 1994; Liu & Rasmusson 1997). We show in this section that the results so far have not conclusively shown advantages of MMs over HH models (except for MMs being able to replicate gating current findings).

Horn & Vandenberg (1984) compared over 20 different MM layouts to find the ideal structure for modelling INa. The model structure identified had a very small likelihood and the difference between the log likelihoods was less than 1 per cent for the top 17 models. The confidence intervals for the parameters of the ‘best fit’ had the same magnitude as the respective parameters, indicating that either all of the proposed models could not replicate the experimental data or the problems in fitting the parameters prevented the success of the method, indicating the importance of parameter identifiability and estimation. Yet, the proposed model was used with only minor adjustments until Clancy & Rudy (2002) proposed a new formulation that was very similar to the original HH formulation with an additional two states for the late or persistent current component.

Other publications discussing the optimal structures of MMs and their differences from HH models used the transient outward potassium current Ito as an example. Campbell et al. (1993b) presented two models, one with an HH structure and one of MM type, and found the dynamics of the models to be equivalent. Liu & Rasmusson (1997) were able to show that the models showed different outputs when an additional state was added to simulate drug block. Ad hoc adjustments of the HH model were not able to reproduce the effects observed in the MM. The conclusion was drawn that HH models and MMs present, in general, different properties when applied to drug block (Liu & Rasmusson 1997).

The open probabilities of the two Ito models are shown in figure 2a,b for a voltage step from −80 to +30 mV. We find small differences that might not be discernable using experimental results due to their inherent noise and variability. But, importantly, the steady-state I–V curves (figure 2c,d) are significantly different as the HH model shows a peak at approximately 10 mV, whereas the steady-state current in the MM increases with increasing voltage. The observed differences are similar to the results found by Liu & Rasmusson (1997) on drug block, which, therefore, can be explained as differences in steady-state current. This means that the two investigated models are different. The additional blocking state, therefore, does not elicit any inherent differences between HH models and MMs. A comparison of the two models using standardized protocols could have elicited the model differences. We will provide such a protocol in §7.

Comparison of two Ito model types. (a,c) HH-type models and (b,d) MMs have been found to be similar in their normal electrophysiological behaviour (Campbell et al. 1993b), but different in their reaction to blocking by a drug (Liu & Rasmusson 1997). On a closer look at (a,b) the open probabilities for a step to +30 mV reveals differences in the dynamics. Even more apparently, (c,d) the steady-state solutions do not match, but show exactly the same differences as observed by Liu & Rasmusson (1997) when adding a blocking state.

4. Method for investigating the identifiability of parameters

The importance of parameter estimation in the modelling process has been discussed in §3. We show in this section that parameters can only be estimated well and uniquely if they are identifiable, and we provide an algorithm to test this in any model.

The algorithm to numerically determine the local identifiability of a set of parameters has been presented and applied previously (Fink et al. 2008a). It is based on results from Reid (1977) on the local information matrix and Cobelli & DiStephano (1980), who generalized this theory and defined the ‘sensitivity function matrix’ (given below).

(a) Sensitivities of the output with respect to parameters

The sensitivity Si,j(t) over time t of a state or an output yi of the system (model) under investigation with respect to a parameter pj is the respective derivative, normalized by multiplying by the parameter and dividing by the state variable—yielding a non-dimensional vector,(4.1)Owing to its dimensionless form, this measure can be used to compare the influence of various parameters (and their uncertainties) on an output. For example, if the sensitivity at a point in time is equal to 1, then an increase in the parameter by 5 per cent would yield an increase in the output by approximately 5 per cent. In the case of the sensitivity being negative, an increase in the parameter would lead to a decrease in the output.

(b) Parameter identifiability using sensitivities

If the output of the system shows the same sensitivities for two parameters, then an increase in one and a decrease in the other would yield the same output as before; this means that the parameters cannot both be identified.

A generalized concept of this is to put all the sensitivities regarding one parameter in a column and compare these with the columns for the other parameters. This leads to the sensitivity function matrix(4.2)where t is the time; si,j, i=1, …, m, j=1, …, k denotes the sensitivities of the system, with respect to the m outputs (measurements) and the k parameters. Assuming that all outputs have the same number of time points n, the resulting matrix has mn rows and k columns.

As explained above, the model parameters are identifiable if and only if the rank of this matrix is equal to the number of parameters to be estimated, or equivalently if and only if the matrix STS is non-singular.

As it is difficult or impossible (dependent on the model) to evaluate the rank of this matrix analytically, we apply a numerical approach, i.e. the matrix is evaluated at fixed nominal parameter values. This yields a criterion for local identifiability, which is a necessary condition for the parameters to be globally identifiable. The software used and further details on sensitivity analysis can be found in the electronic supplementary material.

(c) Numerical methods for solving ordinary differential equations

Unless stated otherwise, all simulations are run using the Matlab (The MathWorks, Natick, MA, USA) solver ode15s, which is a variable-order solver based on the numerical differentiation formulae. Optionally, ode15s could use the backward differentiation formulae (also known as Gear's method), but this is less efficient. Maximum step size was set to 10 ms, because this is the minimum duration of the voltage steps used.

(d) Steady state as initial condition

Voltage-clamp experiments start at a holding potential, i.e. a specific voltage is maintained for a long time to reach steady state. To reduce the computation time for parameter identifiability and estimation, it is advantageous to derive such a steady state analytically.

For the HH model, the steady state for each gate is given or can be easily computed using the transition rates. We developed an algorithm for MMs to derive the steady state for any MM. To be in a steady state, the number of channels changing their confirmation from, for instance, state 1 to state 2 must be equal to the number of channels changing from state 2 to state 1. Given α1 (the probability of a channel in state 1 of changing to state 2) and β1 (the probability of a channel in state 2 of changing to state 1), i.e. , the following must hold for the steady state: S1α1=S2β1 or S2=S1(α1/β1)=S1γ1, where S1 and S2 denote the proportion of channels in states 1 and 2, respectively. This leads to a simple algorithm for determining the steady state by assigning relative proportions (si) to each state according to the equation above. Starting with the first state as s1=1, the relative proportion of state occupancies can be derived by s2=s1γ1, s3=s2γ2, etc. To calculate the absolute probability for each state, one has to divide by the sum of relative proportions (equivalent to the sum of probabilities is equal to one): , where i=1, 2, …, n, and n denotes the total number of states. This algorithm holds if the microscopic reversibility is obeyed (no energy input or output at any of the transitions). An example is given in the electronic supplementary material.

5. In silico voltage clamp protocols

In order to investigate the influence of the data availability on parameter identifiability, three different voltage clamp protocols are used for the in silico experiments (figure 3). Protocol 1 represents the normal range of voltage-clamp steps used to characterize ion channel dynamics (e.g. Berecki et al. 2005). A second protocol (protocol 2) is designed to have the same information content as protocol 1 (see §7), but to be much shorter—it is therefore optimized, but we currently cannot prove that it is the optimal protocol. This protocol optimization is done using generalized sensitivities (Thomaseth & Cobelli 1999). Protocol 3 is very short and (as a proof of principle) does, in general, not provide enough information for parameter estimation. The following voltage clamps form our protocols.

Voltage clamp protocols used to investigate the identifiability of the model parameters. (a) The normal range of electrophysiological voltage clamp experiments to investigate activation, deactivation, inactivation and reactivation (total length 39.8 s). (b) A shortened protocol that contains all the information of the long protocol shown in (a) (13.5 versus 39.8 s, respectively). (c) The short (2.8 s) voltage clamp protocol does not provide sufficient information to identify as many parameters as (a, b).

(a) Protocol 1

Extensive protocol (figure 3a); (i) activation and inactivation: repeatedly from a holding potential of −80 mV (for 600 ms) to one voltage from −60 to +40 mV (with 20 mV increments) for 1000 ms and back to the holding potential; (ii) deactivation: repeatedly from a holding potential of +60 mV (for 1000 ms) to one voltage from −100 to 40 mV (20 mV increments) for 1000 ms and back to the holding potential; (iii) reactivation: repeatedly from a holding potential of +60 mV (for 1000 ms) to one voltage from −100 to +20 mV (20 mV increments) for 10 ms, then up to +40 mV for 100 ms and down to −60 mV for 200 ms and back to the holding potential.

(b) Protocol 2

Shortened protocol (figure 3b); (i) activation and inactivation: from a holding potential of −80 mV (600 ms) to −60, −40 or +20 mV for 1000 ms and back to −80 mV; (ii) deactivation: from a holding potential of +60 mV (1000 ms) to −100, −70 or −40 mV for 1000 ms and back to holding potential; (iii) reactivation: from a holding potential of +60 mV (1000 ms), a voltage step to −100 or −80 mV for 10 ms, then up to +40 mV for 100 ms and down to −60 mV for 200 ms and back to the holding potential.

(c) Protocol 3

Very short protocol (figure 3c); from a holding potential of −80 mV, a voltage step to +30 mV for 800 ms, step down to −120 mV for 600 ms, step to 0 mV for 20 ms, step to −100 mV for 400 ms, to −80 mV for 300 ms, to −20 mV for 300 ms and back to −80 mV for 300 ms.

6. Parameter identifiability in voltage-gated ion channel models

Theoretical results on single channel experiments have shown that MMs can theoretically be identified uniquely using a single channel voltage step clamp (Fredkin et al. 1985; Bauer et al. 1987). However, with real, noisy data, the reliable analysis of first latencies is almost impossible (Magleby 1992). Therefore, we investigate the identifiability of parameters based on whole-cell measurements of ion channel kinetics.

This section discusses the necessary conditions for obtaining a unique parameter set that best fits the available experimental data, assuming that a model structure is provided with guesstimates for the parameters. We examine a number of voltage-gated ion channel models (see §2) to investigate the identifiability of their parameters from whole-cell currents using electrophysiological voltage clamp experiments (see §5). Unidentifiable parameters make it difficult to use numerical algorithms on a computer to identify the optimal parameters, because different initial parameter guesses will lead to different results and the search scheme will be very slow due to numerical problems such as inversion of the (close to) singular Jacobian matrix and the number of local minima.

There exist two kinds of non-identifiable parameters: a priori and numerically unidentifiable. A priori unidentifiable parameters can be found directly from the model equations. They are independent from the actual parameter values and the given experimental data. The most common example found in the models investigated is given by the following equation, where the parameters p1, p2 cannot both be identified:(6.1)Any parameters and yield the same solution for this equation with any parameter c. Thus, there exists no unique set of parameters. This issue can be resolved by either setting p1=1 or p2=0. Fixing one of the parameters during the parameter estimation process would work as well, but would be misleading in two ways: (i) it would indicate that there is further information available on that parameter (if there actually were information, one would not have to estimate this parameter) and (ii) people reusing the model structure would not be aware of the fact that this one parameter is not estimated numerically.

Out of 13 models discussed in this work, more than half of the models exhibit a priori unidentifiable parameters (see table 1 and the electronic supplementary material for specific parameters).

Parameter identifiability using standard voltage clamp experiments. (Comparison of published ion channel models. Unidentifiable parameters are given as a priori or ‘numerically’. A priori identifiability is dependent exclusively on the model formulation. The number of numerically identifiable parameters is dependent on the available experimental data (compare protocols) and its quality (see text for details). Note that protocols 1 and 2 provide similar parameter identifiability.)

‘Numerical unidentifiability’ is much harder to define generally. Parameters may or may not be identifiable depending on the size and quality of the available experimental data and distance of the initial parameter guess to the optimal (but unknown) parameters of the model. There is also a large influence of the experimental protocol on the ‘numerical identifiability’ of the model parameters, e.g. two measurements at steady state allow the estimation of only two parameters. Therefore, we investigate three different experimental voltage-clamp step protocols (see §5 and figure 3; assuming no or little noise in the output): a long one, consisting of a complete set of experiments to investigate activation, deactivation, reactivation and inactivation (protocol 1); a shortened one (protocol 2, see §7 for details); and a very short one (protocol 3).

The results given in table 1 show that in only 4 of the 13 models investigated the parameters may be identified uniquely. However, this does not imply that in these four models the parameters actually can be identified at a low computation cost and that there are no local minima. At least the model parameters have been shown to be locally identifiable, which is a necessary prerequisite for global identifiability. Considering that the a priori unidentifiability can be fixed as mentioned above, we could add adjusted formulations of the IKr models by Mazhari et al. (2001) and ten Tusscher et al. (2004) to the set of locally identifiable models.

The analysis provides further insights into which parameters are unidentifiable (see tables in the electronic supplementary material), such that either these parameters can be neglected (set to zero) or the equations can be changed to resolve the problem. For instance, for the INa model by Clancy & Rudy (2002), the inactivation rate constant is given as(6.2)Even in the case of large voltage values, the second part of equation (6.2) is very small compared with the first part, and therefore the respective parameter is flagged when running our sensitivity analysis software as ‘numerically unidentifiable’.

The threshold for numerical identifiability cannot be well defined; therefore, the actual number of unidentifiable parameters can vary. There are cases where the identifiability problem can be easily identified—as the one given in equation (6.2). Another possibility is that a gap appears in the calculated eigenvalues, indicating that the parameters after the gap cannot be identified (Demmel & Kagström 1993). If there is no obvious gap, but the eigenvalues span a wide range (bad condition of the sensitivity matrix), then the model has too many parameters and the structure has to be changed in order to reduce the number of parameters.

Note that the three published INa models we consider are all overparametrized. In order to show that one can construct a model for INa that is identifiable, we use the structure/layout of the model by Clancy & Rudy (2002; figure 1b), exchange the opening and closing transitions with HH-type transitions (thus removing seven parameters) and remove the two a priori unidentifiable parameters in these transition rate formulations to obtain(6.3)(6.4)(6.5)(6.6)We also remove the numerically unidentifiable parameters from the equation of β3 (see equation (6.2) described above) and leave the other equations unchanged,(6.7)(6.8)(6.9)(6.10)The identifiability analysis in table 1 indicates that for this model (FinkINa) the remaining 15 parameters could all be identified using any of the three presented protocols.

For the ICaL model by Faber et al. (2007), the algorithm identifies one numerically unidentifiable parameter. However, they also set ωfs=Φs without further explanation (i.e. the transition rate from the open state and the fast inactivated state to the slow inactivated state are set to be equal). Without this assumption, two more parameters have to be estimated.

The parameters of the two Ito models by Campbell et al. (1993a) are identifiable, which are coherent with their rather small number of parameters (6 and 9, respectively), whereas the model by Greenstein et al. (2000) has unidentifiable parameters. There is also one numerically unidentifiable parameter in the IKs model by Silva & Rudy (2005).

We are not aware of a publication that has investigated the question of necessary or sufficient data for estimating the parameters of voltage-gated ion channels from whole-cell currents.

Most models are fit with as much data as possible using a weighted cost function for the error minimization (e.g. Greenstein et al. 2000; Fink et al. 2008b), i.e. they are ‘developed on the basis of a wide range of data presented in the experimental literature’ (Clancy & Rudy 2002). This means that not all data will have been collected using the same conditions (such as bath and pipette solutions, etc.), and thus may not be compatible. Using all data also means that there are no additional data left for model validation. Furthermore, additional uncertainty for parameter estimates is added when redundant information is added to the data (Banks et al. 2007).

We compare three different voltage-clamp step protocols (as described in §5) to investigate the information content of protocols of different lengths (2.9, 13.5 and 39.8 s for protocols 1, 2 and 3, respectively).

Protocol 1 (figure 3a) denotes the normal set of experimental protocols used to investigate voltage-gated ion channel kinetics. Adding more voltage-clamp steps (e.g. to investigate the reactivation of the current) did not change any of the presented results (not shown).

The second protocol (protocol 2) is designed using generalized sensitivities (Thomaseth & Cobelli 1999). Figure 4a shows the generalized sensitivities of the IKr model by Wang, Q. et al. (1997) for the extensive physiological protocol 1. Flat lines or repetitive patterns indicate that there is no additional information content in these parts of the protocol. By removing these parts from the protocol, one obtains the shortened protocol 2. The generalized sensitivities of protocol 2 do not show any flat lines or repetitions anymore (figure 4b), so it is optimized in that respect.

The generalized sensitivities denote the increase in information for each time point. Horizontal lines indicate no substantial information gain, whereas steep curves denote important parts of the protocol for the respective parameter. The generalized sensitivities for four parameters of an IKr model for (a) protocol 1 (figure 3a), which shows some length of horizontal lines and very similar steps at the end, and (b) a shortened protocol (figure 3b) without horizontal or repetitive parts.

To see if protocol 2 is still equivalent to protocol 1, we use the methods of parameter identifiability described in §§4 and 6, as they allow us to determine which experimental protocols contain enough information to estimate the same amount of parameters as the long experimental voltage step protocol (protocol 1).

The results on numerical identifiability (see table 1) show that the shortest protocol (3) cannot identify as many parameters as protocol 1 or 2. There are only two models that could be fitted using this protocol (FinkINa and the Ito MM by Campbell et al. 1993a). Therefore, it is not suitable as a general protocol.

Protocols 1 and 2 contain similar information, as can be seen from the same number of identifiable or unidentifiable parameters for these protocols in table 1; more importantly protocol 2 is much shorter in length. With more rigorous methods, it might be possible to reduce this protocol further, but it is already very comprehensive and the results (table 1) show that it can be used to fit a wide range of models for voltage-gated ion channels, such as INa and (the voltage-gated part of) ICaL, Ito, IKs and IKr.

8. Computational effort to solve ion channel model equations

The computational effort to solve ion channel model equations varies a lot between models. The parameter identification analysis for the model by Clancy & Rudy (2002) took approximately 7.8/3/1.5 min for protocols 1/2/3, respectively (on an AMD Turion 2.0 GHz), whereas the analysis for the model by Irvine et al. (1999) took approximately 96/36/15 hours, respectively. This remarkable difference in computation time between similar models underlines the importance of computational speed. Similar results can be expected for parameter estimation and also for multicellular whole ventricle simulations.

In this section, we investigate the difference in importance between the number of states and the stiffness of the models, which clarifies the appropriateness of MMs for multiscale modelling. Then, we discuss ways to improve the computational speed when considering voltage-clamp step protocols, for instance for parameter estimation.

(a) Computational effort to solve ion channel models

To investigate the computational efficiency of ordinary differential equation (ODE) models, two solution methods for simulating voltage clamp protocol are used and compared (table 1) using voltage clamp protocol 2 as the test case. Specifically, the solutions are obtained (i) by restarting the simulation at each change in voltage and (ii) solving over the whole time interval.

For solutions using an ODE solver, the computational cost is dependent on (i) the cost for deriving the approximation of the solution to an ODE at any given time point and (ii) the actual step size taken by the variable step size solver.

The first is related to the computational cost at each time point, which is dependent on the efficiency of calculating the derivatives of the states, and the Jacobian of the system. The second is associated with the stiffness of the system, which correlates with the large differences between the time constants of the various state variables. (Physiologically, this is related to, for instance, the rapid inactivation of INa and its comparatively slow reactivation process. The speed of these two processes can be orders of magnitudes apart.)

To investigate whether the number of states is an important determinant for the computation time, we compared the results for the HH model by Priebe & Beuckelmann (1998) with its equivalent MM formulation (denoted with * in table 2). In the case of Priebe & Beuckelmann's (1998) model, the system with three HH-gating variables was written in the MM formulation using all 16 states. Such a representation is equivalent since the three m, h and j gating variables in the HH formulation of INa represent, respectively, 4, 2 and 2 states, which when coupled produce 4×2×2=16 states in the MM case. The systems are analytically equivalent, thus the stiffness has to be the same. Therefore, the difference in the simulation times of 2.27 s for the HH model compared with 4.00 s for the MM comes from the increase in the number of states.

Speed, stiffness and number of states. (Comparison of published and modified (asterisks) models. The modified Priebe model uses the MM formulation instead of the HH one. The modified Faber model uses the HH model for the calcium-dependent inactivation. The table columns present model, respective current, type of current formulation, number of state variables, number of actual time constants, ODE solver time (where the first value corresponds to a piecewise solution and the second value corresponds to a normal solution) and exponential solver time. The simulation time is determined by both the number of states and the stiffness of the equations. The stiffness determines the size of the variable time steps (correlating with ODE solver time), whereas the size of the system determines the calculation time of the Jacobian at each time step (correlating with the time for exponential solutions).)

As a second example, we decreased the number of states of the ICaL model by Faber et al. (2007) by using a model with mixed formulation of HH model and MM (denoted with asterisks in the table). In this model, the calcium-dependent inactivation of the model is completely independent of all other states, and can therefore be modelled using an HH formulation. Thereby, the number of states can be reduced from 14 to 8, decreasing the simulation time from 3.35 to 2.86 s. The resulting model can be described as follows (see details in the electronic supplementary material):(8.1)with(8.2)where OICaL, OV and OCai denote the open probabilities of the channel, the voltage-dependent gating and the calcium-dependent inactivation, respectively. δ and θ are the transition rates for inactivation and reactivation of the calcium-dependent inactivation as given in Faber et al. (2007).

However, in both cases, the computation time of the model with more states is less than double that of the smaller model.

The second parameter influencing the computation time is the stiffness of the system. This can be investigated comparing the results for the ODE solver with and without resetting at each voltage step (first/second number in table 2, column 6). The investigated models show that MMs can go from a non-stiff state (Silva & Rudy 2005) to a very stiff one (Irvine et al. 1999) as can be seen in the simulation times for the results using an ODE solver 1.42/1.89–2.98/16.4 s, respectively. This is already a fivefold difference, which can be enhanced when more complex simulations are done, similar to our sensitivity analysis, where we get a 600-fold difference between the INa models by Irvine et al. (1999) and Clancy & Rudy (2002).

(b) Reducing the computation time by using analytical solutions

The computation times in table 2 look rather short, but as the solution for the voltage-clamp protocol has to be calculated over and over again during parameter estimation, a sevenfold speed-up means the reduction of computation time from a week to a day.

As the voltage is piecewise constant for these protocols, one can use analytical solutions for each voltage step (see the electronic supplementary material for details) instead of an ODE solver to speed up the computations. The last column in table 2 shows the computation times using the analytical solutions. The speed-up is of the order of 3–10 times for MMs, approximately 16 times for the mixed (MM and HH) ICaL model, but more than 40 times for the HH model types.

Note that using the same model, but using different formulations and solution algorithms, can lead to computation times of 4.0 versus 0.038 s, which is a 100-fold difference (compare last column for Priebe98 and first value of sixth column for Priebe98* in table 2). However, this significant difference holds only for the HH model, due to the costs of deriving the eigenvalues and eigenvectors of the MM matrix.

9. Conclusions

The goal of mathematical modelling in electrophysiology is to develop mathematical descriptions of ion channels and whole cells to provide further insights into underlying mechanisms and to predict the outcome of experimental results in order to inform future research.

(a) Versatility in model structure

In this paper, we focus on the modelling of voltage-gated ion channel kinetics. The derivation of MMs for exchangers and pumps can more easily be based on biophysical properties and binding kinetics (see Smith & Crampin 2004), and thus is more robust, while determining the structure of voltage-gated ion channel kinetics is challenging as demonstrated in this paper (§3). This raises the question as to whether the model is able to provide insights into the molecular structure of the ion channel. We show in §3 that past efforts to determine the optimal model structure (e.g. in Horn & Vandenberg 1984) have failed owing to problems with parameter estimation and proper model comparison.

Furthermore, currently, MM rates often contain only linear terms for the voltage-dependent Gibbs free energy G(V). It has been shown that the quality of fit can be improved without the need of additional states by considering formulations that include, for instance, (i) nonlinear terms (Destexhe & Huguenard 2000) or (ii) multiple energy barriers (Ozer 2007), thus providing even more options to consider for an ‘optimal’ model.

(b) Parameter identifiability of voltage-gated ion channel models

In order to obtain a more robust parameter estimating approach and the best layout of the ionic model, we investigate the parameter identifiability and analyse different experimental voltage-clamp step protocols.

Whereas one could simplify the estimation of parameters by using different protocols to fit different parameters (Vecchietti et al. 2006), it is critical to fit all of the model parameters together as has been shown by various authors (Willms 2002; Lee et al. 2006).

We find that 9 out of 13 ion channel models have unidentifiable parameters (see §6) that make unique parameter estimation impossible. Moreover, all of the three considered published INa models are overparametrized, so we provide an improved model formulation with identifiable parameters.

(c) Newly developed short, but sufficient voltage-clamp step protocol

The shortened voltage-clamp step protocol presented in §7 has the same information content as the extensive protocols normally used to investigate voltage-gated ion channel kinetics. This protocol can therefore be used as a fingerprint for the ion channels investigated (INa, ICaL, Ito, IKs and IKr). Furthermore, this is the first time that the question has been investigated as to how much experimental data are necessary for fitting a voltage-gated ion channel model from whole-cell currents; the remaining data can be used for model validation.

(d) Computational speed: stiffness versus number of states

The structure and especially the number of states in MMs have often been thought to slow down the computational speed of multicellular simulations. We show (see §8) that the decrease of the maximum step size due to stiffness affects the computational costs significantly more than an increase in the number of states.

However, the stiffness of a model can be reduced by using a quasi-steady-state approach for the fast components in the model, assuming that the very fast transitions are always at equilibrium (see Hinch et al. 2004; Smith & Crampin 2004).

The analytical solution of voltage-clamp step protocols decreases computation time by up to 90 per cent for MMs (by 97.5% for HH models) compared with ODE solutions and should therefore be used in parameter estimation.

(e) Cell-to-cell variability and differences between species

Current models show a huge variety in dynamic and steady-state behaviour as shown in figure 5 for the INa and the IKr currents. Such a variability is based more on the parameter identifiability problems specified in this paper rather than on the real cell-to-cell variations or species differences. Results using the optimized experimental voltage clamp protocol presented can be obtained from individual cells, yielding information on cell-to-cell variabilities and covariances of parameters (this has been done already for IK1 in Fink et al. (2008b)).

Comparison of model results for the optimized voltage clamp protocol. Apparently, there is a large variation between models for INa in (a) as well as for IKr in (b). This variation is partly due to species differences and variations in experimental data. To investigate how much of this is due to cell-to-cell variability, the models (HH or Markov) have to be fit to the results from just one cell, i.e. all parameters have to be identifiable from the experimental output of one cell, for instance, by using our shortened voltage clamp protocol (figure 3b). This protocol also provides a unique fingerprint for each voltage-gated ion channel.

The set of tools we use in this paper provide the possibility for higher standards in model development and can help to obtain automated methods for model development, especially to identify structure for each ion channel and to further compare advantages of HH and Markov modelling approaches.

Acknowledgments

The work in Oxford is supported by grants from the EU BioSim Consortium (No. 005137), the EU preDiCT grant (ICT-2007-2-224381) and the Novartis Pharma AG.

Footnotes

One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications II’.