About this Journal Submit a Manuscript Table of Contents
Computational and Mathematical Methods in Medicine
Volume 2013 (2013), Article ID 706195, 20 pages
http://dx.doi.org/10.1155/2013/706195
Research Article

Optimisation of a Generic Ionic Model of Cardiac Myocyte Electrical Activity

Graduate School of Biomedical Engineering, University of New South Wales, Sydney, NSW 2052, Australia

Received 19 November 2012; Revised 11 February 2013; Accepted 18 February 2013

Academic Editor: Henggui Zhang

Copyright © 2013 Tianruo Guo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A generic cardiomyocyte ionic model, whose complexity lies between a simple phenomenological formulation and a biophysically detailed ionic membrane current description, is presented. The model provides a user-defined number of ionic currents, employing two-gate Hodgkin-Huxley type kinetics. Its generic nature allows accurate reconstruction of action potential waveforms recorded experimentally from a range of cardiac myocytes. Using a multiobjective optimisation approach, the generic ionic model was optimised to accurately reproduce multiple action potential waveforms recorded from central and peripheral sinoatrial nodes and right atrial and left atrial myocytes from rabbit cardiac tissue preparations, under different electrical stimulus protocols and pharmacological conditions. When fitted simultaneously to multiple datasets, the time course of several physiologically realistic ionic currents could be reconstructed. Model behaviours tend to be well identified when extra experimental information is incorporated into the optimisation.

1. Introduction

The electrical activity of cardiac myocytes, including action potential (AP) waveforms and underlying membrane currents, has been extensively studied using a combination of microelectrode recording and mathematical modelling techniques [1]. The latter in particular has been utilised to provide a more quantitative and integrative understanding of the underlying mechanisms of cardiac electric activity. Biophysically detailed ionic models of cardiac cell electrophysiology are able to accurately reproduce a large range of behaviours, including membrane potential waveforms, specific ionic currents under voltage-clamp protocols, AP membrane current dynamics, and Ca2+ alternans [2, 3]. However, the ever-increasing number of variables in such models renders them computationally expensive when integrated into higher-dimensional tissue or whole heart simulations, limiting their utility.

As an alternative, simplified phenomenological models have been widely utilised in electrophysiological simulations due to their minimal complexity and computational cost. The first such generic model was the Fitzhugh-Nagumo (FHN) formulation published in 1961 [4], simplifying the four-variable Hodgkin-Huxley (HH) nerve axon model [5] into a two-variable formulation by eliminating gating variables with rapid time constants. Since its publication, FHN-type models have been frequently used in multicellular tissue simulations [69]. In 1998, Fenton and Karma published an improved phenomenological model based on the biophysically detailed Luo and Rudy [10] and Beeler and Reuter [11] ventricular cell ionic models in order to simulate ventricular fibrillation [12]. Given appropriate parameters, the ability of the three-variable Fenton-Karma model to reproduce AP duration restitution properties was shown to be comparable to biophysically detailed models [13]. By introducing an additional variable to this model, the Fenton-Cherry modification was also able to reproduce AP waveforms from pulmonary vein and left atrial myocytes [14]. Despite the successful application of these phenomenological models, their inherent oversimplicity may restrict their utility, particularly when simulating complex phenomena such as electrical remodelling during sustained arrhythmia [15] or the effects of selective ionic channel blockers [16]. For example, it is doubtful whether such models can accurately reproduce the range of AP waveforms recorded from the same myocyte under various electrical pacing frequencies or multiple degrees of selective ion channel block by drugs, due to the small number of membrane currents in these models.

Another significant challenge in cardiac single cell ionic modelling lies in estimating the host of parameters governing the kinetics and densities of the various membrane ion channels, pumps, and exchangers, as well as parameters governing intracellular ionic cycling and buffering mechanisms, all from readily available experimental data obtained in a given myocyte. Given enough parameters, ionic models may only be able to accurately reproduce ionic mechanisms under the precise experimental conditions they were fitted to in the first place. Models with modified parameters to fit another set of experimental data may lose their original mechanisms [3, 17], limiting their predictive utility. Moreover, the ever-increasing complexity of such models makes parameter estimation a highly difficult and time-consuming task. Although previous studies have undertaken parameter optimisation in cardiac ionic models [18, 19], these have all used a relatively limited subset of parameters (typically maximum membrane conductances only) to fit AP records. Until now, a very limited number of optimisation algorithms have been successfully used for large-scale optimisation of cardiac cell ionic models.

In response to these challenges, we have formulated a simplified HH-type generic model to accurately reconstruct a range of AP waveforms recorded from tissue-intact myocytes in rabbit sinoatrial (SAN) and right and left atrial (RA, LA) tissue preparations. The model structure is flexible and modular and the optimised models are able to reproduce complex behaviour such as the change in AP morphology due to selective ion channel block or high-frequency paced stimulation. Model parameters were estimated using a computationally simple, multiobjective AP optimisation approach based on a custom curvilinear-gradient method [20]. The improvement in membrane current reconstruction gained by optimising the model to simultaneously fit different AP experimental waveforms was examined.

2. Methods

2.1. Generic Cell Model Formulations

The generic form of the model is given by where (mV) denotes the transmembrane potential, denotes time, (μF/cm2) is the specific membrane capacitance, (nA/cm2) refers to the time-independent (leakage) current, (nA/cm2) is the applied stimulus, (nA/cm2) denotes the th time-dependent ionic current density, and is the user-defined number of time-dependent currents, the value of which depends on the complexity of the data to which the model is to be fitted. The ability to set the number of currents makes the model generic and modular. The time-dependent currents are represented by a two-gate HH scheme: where (μS/cm2) is the maximum ion channel conductance, and are dimensionless gating variables, and (mV) is the reversal potential for membrane current . The time-independent leakage current is described by where (μS/cm2) and (mV) are the maximum ion channel conductance and reversal potential, respectively. For the gating variables in (2), HH kinetics [5] are employed to specify first order differential equations (ODEs): where , and (s−1) represent opening and closing rates for the corresponding gating variable. These rates are given by sigmoidal functions of the membrane potential : where (s−1), (mV−1), and (mV) are parameters specific to each rate, with each parameter value shared between the and pair for each gating variable.

Therefore in total, the generic model contains ODEs and parameters. In principle, the higher the value of , the more degrees of freedom and the more complex are the electrophysiological behaviours which can be reproduced, at the cost of decreased computational efficiency.

2.2. Parameter Optimisation

A custom curvilinear gradient-based least-squares optimisation method, combining the advantages of both Newton and steepest-descent methods [20], was carried out on a standard desktop PC using Matlab software (The Mathworks Inc., USA). The curvilinear gradient method may be briefly described as follows: we assume that an array of data points are to be fitted by an ODE system , where is an array of parameters. The residual array between model and data is given by In general, the model output will be a nonlinear function of the parameter array p. However, the residual can be linearized using the approximation: where is the array of parameter increments relative to the current parameter set , and is the array of residuals. A least-squares scalar objective function, , to be minimised is given by where denotes the transpose, is the Jacobian matrix, is the current objective value, denotes the gradient of the objective at the current parameter point , and is the Hessian matrix. If the model is well determined in its parameters, will be nonsingular at the global minimum for .

Equation (8) represents a quadratic form in the parameter increment vector . The exact minimum of this quadratic form will occur at

and represents a full Newton step to the minimum objective value. In practice, however, the least squares scalar objective given by (8) will not be accurate for large , where the linear approximation of the residual in (7) breaks down. When is sufficiently small, (8) will be accurate. We can however search for the minimum objective along a curvilinear trajectory through the -dimensional parameter space given by the curve of steepest descent of the quadratic form (8). This trajectory is given by where is a parameter between 0 and infinity, with corresponding to the current parameter location and corresponding to the full Newton step mentioned previously. When the objective minimum has been found along the trajectory, the quadratic form (8) is recomputed, and another search is performed along a new curvilinear path. The approach is iteratively carried out until the least square objective is locally minimised. Full details of the method are described in Dokos and Lovell [20] including random restarts and iterative reweighting to find the global least squares minimum. In addition, some of the experimental AP waveforms recorded in the present study exhibited stimulus artifacts which could confound the optimisation process. These were effectively removed by assigning a weight of zero to the residual between model and data in these regions.

When optimising the model to fit multiple datasets simultaneously, the residual in (6) was formed by appending together the residuals of the model and the corresponding individual datasets. As a result, the calculation of the Jacobian matrix in (7) was slightly modified. Three cases can be considered.(i)Multiple data ( datasets) fitted using the assumption that each dataset shares the same parameter values, for example, fitting multiple AP data recorded from the same cell in response to different pacing conditions. will be of size , where is the total number of data points across all records and is the number of optimised model parameters.(ii)Multiple data ( datasets) fitted using the assumption that each model uses a unique set of parameters to fit each experimental dataset. Data-specific parameters to , each of size , are used for datasets 1 to . This process is equivalent to performing multiple single dataset optimisations independently, and then will be of size . (iii)Multiple data ( datasets) fitted using a combination of both shared and data-specific parameters. This is the case for the drug-specific and tissue-specific optimisation scenarios presented. For data-specific parameters (i.e., parameters unique to each dataset), will be of size .

Compared with single dataset fitting, more computational resources are required for multiobjective optimisations due to the larger size of the Jacobian matrix, as well as the fact that more local minima are likely to be present in the objective parameter space.

2.3. Cardiac Tissue Recordings

New Zealand White rabbits (6–24 months old) were anesthetised with 5% isofluorane, and 1000 IU of heparin was administered intravenously. A thoracotomy was performed and the heart was rapidly excised and placed in cold cardioplegia solution. SAN-RA or LA appendage tissue preparations were dissected and placed in a recording chamber, superfused with Tyrode’s solution and oxygenated with 95% O2 and 5% CO2 to maintain the pH at ~7.4. Intracellular APs were recorded using sharp glass microelectrodes (resistance 50–100 MΩ). Recordings were amplified (gain × 10) and filtered (low pass, cutoff 10 kHz) using an Axoclamp 2B amplifier (Axon Instruments, USA) and sampled at 20 kHz using a USB-6251 analog/digital converter (National Instruments, USA) and a custom build data acquisition software programmed in Labview (National Instruments, USA). For the LA experiments, the tissue was paced using a STG1002 stimulator (MultiChannel Systems, Germany). Monophasic suprathreshold pulses, 2 ms in duration, were delivered using Teflon coated bipolar stainless steel electrodes (125 μm diameter). All experiments were conducted in accordance with Australian animal research guidelines and were approved by the University of New South Wales Animal Ethics and Care Committee.

3. Results

3.1. Results of Optimised Minimal Generic Model

A minimal generic HH-type ionic model with two time-dependent membrane currents, one inward and one outward, in addition to a background leakage current was fitted to three consecutive spontaneous APs recorded from central sinoatrial node (cSAN), peripheral sinoatrial node (pSAN), and RA tissue-intact myocytes from rabbit sinoatrial tissue preparations. Optimised APs along with the corresponding ionic currents are shown in Figure 1. Lists of optimised parameter values and initial values for all model variables for each cell are given in Tables 1 and 2. In accordance with the data, the pSAN model-generated APs exhibited a faster upstroke, a more positive overshoot, and more negative maximum diastolic potential compared to the cSAN. Both types of SAN APs exhibited a slow depolarisation (pacemaker) phase which was absent in the RA cell, and, as a result, stimulus pulses of 2 ms duration and 30 μA/cm2 amplitude were applied to trigger APs in the RA model. There was a gradual transition in the shape of the AP and ionic current waveforms from central SAN to atrial tissue. In particular, the transient spike of the inward current (ionic current 1) was nearly absent in the cSAN model and progressively increased both in magnitude and upstroke rate from pSAN to RA. The root mean square (RMS) error between the optimised models and data was 1.02 mV, 1.53 mV, and 1.69 mV for cSAN, pSAN, and RA myocytes, respectively.

tab1
Table 1: Parameter values for single AP dataset-based optimisation.
tab2
Table 2: Initial variable values for single AP dataset-based optimisation.
fig1
Figure 1: Optimisation based on single AP dataset. (a) Top panels: optimised model generated APs overlaid with experimental APs recorded from central sinoatrial node (cSAN), peripheral sinoatrial node (pSAN), and right atrial (RA) intact myocytes from rabbit sinoatrial tissue preparations. Lower panels: corresponding ionic and leakage currents generated by each cell model. (b) Phase plot ( versus ) of model APs overlaid with experimental data. The experimental voltage derivatives were obtained from first order finite differencing of the membrane voltage data.

To assess the generalizability of the HH generic model and optimisation procedure in fitting a wide range of AP waveforms, five separate spontaneous AP recordings from each of the previous myocyte types were fitted using the three-current version of the generic model. Stimulus pulses of 2 ms duration and variable amplitude were used to excite the atrial models. The five sets of data comprised one group described previously (Group 1), plus an additional four sets (Groups 2–5) for each of the three myocyte types (cSAN, pSAN, and RA). Despite the inherent variation between APs recorded from the same myocyte type in different tissue preparations, the generic model was able to fit each dataset using only an inward, an outward, and one background membrane current (Figure 2). The RMS error for each cell type (, mean ± standard deviation) was  mV (cSAN),  mV (pSAN), and  mV (RA). Compared with the fits shown in Figure 1 (Group 1), the mean RMS is marginally increased for all three cell types, likely due to the fact that unsmoothed experimental datasets (with mean peak-to-peak noise levels of ±0.96 mV) were used in fitting data Groups 3–5. Nonetheless, the model was able to reproduce the variability in AP waveform from the same cell type in different preparations.

706195.fig.002
Figure 2: Generic model fits to a range of APs recorded from many myocytes. Optimised models overlaid with five recorded intact tissue APs from each of rabbit central SAN (cSAN), peripheral SAN (pSAN), and right atrial (RA) myocytes. Each single dataset was separately fitted with two time-dependent membrane currents and one leakage current. The average and standard deviation of the root mean square (RMS) error of the fits is also shown.
3.2. Multiobjective Action Potential Optimisation
3.2.1. Multidataset Optimisation with Shared Parameters (Uniformly Paced Left Atrial Data)

The generic model was also optimised to simultaneously fit two morphologically different LA AP waveforms from the same cell, recorded in response to stimulation at pacing intervals (PIs) of 400 ms and 200 ms (Figure 3(a)), using a single set of parameters. A total of five time-dependent ion currents and one leakage current were required to simultaneously fit the two sets of data, with results shown in Figure 3. Stimulus pulses of 2 ms duration and 13 μA/cm2 in amplitude were used to elicit APs in the model. For a PI of 200 ms, the model generated AP characteristics and corresponding time-dependent ionic currents revealed beat-to-beat variations. Parameter values and initial values of all model variables are listed in Tables 3 and 4. The RMS errors between the optimised model and corresponding experimental data were 2.01 mV (PI = 400 ms) and 3.22 mV (PI = 200 ms). This optimised model was then used to predict the APs elicited at a PI of 300 ms, a dataset which was not used in the optimisation process. The additional experimental dataset can be reproduced with an RMS error of 2.46 mV (Figure 3, PI = 300 ms).

fig3
Figure 3: Multiple-dataset optimisation and validation of generic model to fit APs in response to uniform pacing at different frequencies. Experimental APs were obtained by pacing a left atrial intact myocyte at three different pacing intervals (PIs). Model optimisation was repeated three times (runs 1, 2, and 3), each starting at randomized initial parameter values. AP fits obtained for each run were very similar. (a) Three groups of optimised AP fits in response to pacing at intervals (PIs) of 400, 200, and 300 ms. The generic model with five time-dependent currents and one leakage current was simultaneously fitted to the first two datasets (PI = 400 ms and 200 ms) using a single set of model parameters. The optimised model was validated by its ability to reproduce AP responses to pacing at a PI of 300 ms, a dataset not used in the model optimisation. (b) Plots of model generated time-dependent currents for each pacing protocol for each optimisation run. Note the marked AP and membrane current beat-to-beat variations at PI = 200 ms.

To test the reliability of the previous multiple-dataset based optimisation, two additional optimisation runs were carried out utilising the same LA paced data mentioned previously. Both optimisation runs started with the same stimulation settings, but used randomly generated initial parameter values. Iterations were terminated when the RMS error between model and experimental data reached a similar value to that obtained earlier. The resulting RMS errors between the optimised models and corresponding data were 2.01 and 2.40 mV (PI = 400 ms), 3.25 and 3.26 mV (PI = 200 ms), and 2.58 and 2.91 mV (PI = 300 ms) for optimisation runs 2 and 3, respectively.

Although model-generated APs from each optimisation run are almost identical, parameter values displayed significant variation between each run (see Tables 3 and 4 for a list of parameters and initial values for model parameters). The maximum relative difference (as a percentage of mean parameter value) is approximately 170% and the mean relative difference is 19%, indicating that model parameters cannot be uniquely identified, even when they are obtained from multiple dataset optimisations. Interestingly, despite the randomness of the initial parameter values for each optimisation run, the corresponding time-dependent ionic current waveforms are nearly identical and physiologically reasonable, except for differences in their scaling (Figure 3(b)).

tab3
Table 3: Parameter values for LA multidataset-based optimisation.
tab4
Table 4: Initial model variable values for LA multidataset-based optimisation.
3.2.2. Multidataset Optimisation with Shared Parameters (Randomly Paced Left Atrial Data)

A rabbit left atrial (LA) tissue preparation was electrically paced at uniform frequencies with a pacing interval (PI) of 200 ms and 400 ms, as well as randomly paced, and APs responses were recorded from the same myocyte for each pacing protocol. The pulse amplitude and duration were fixed for all protocols. A sequence of 100 random pulses was generated from a normal distribution of PIs, with mean and standard deviation of 275 ms and 69 ms, respectively. A sequence of seven pulses was selected and used in optimisation. In addition to AP alternans observed at a PI of 200 ms, more significant beat-beat variations in AP morphology were demonstrated in the randomly paced dataset, as shown in Figure 4 (lower panel).

fig4
Figure 4: Multiple-dataset generic model fits to APs recorded in response to uniform and random pacing protocols. From top to bottom: AP fits in response to uniform pacing at intervals of 400 and 200 ms and a sequence of random pacing intervals (PIs). The random sequence was generated from a normal distribution of mean 275 ms and standard deviation 69 ms. The generic model with seven time-dependent currents and one leakage current was simultaneously fitted to all three experimental datasets using a single set of model parameters.

Using the generic ionic model, a total of seven time-dependent membrane currents and one leakage current were required to fit the multiple datasets simultaneously using a single set of parameter values (Table 5 lists the shared parameters), with distinct variable initial values for each dataset (Table 6). Stimulus pulses of 2 ms duration and 30 μA/cm2 in amplitude were used to excite the model at the appropriate PI. The optimised model was able to accurately reproduce the experimental AP waveforms, even for the random-paced protocol (Figure 4). RMS errors between the optimised model and corresponding experimental data were 2.94 mV (PI = 400 ms), 3.51 mV (PI = 200 ms), and 3.91 mV (randomly paced). Figure 5 illustrates the corresponding membrane currents generated by the optimised model when paced with the three protocols. Note that the generic model structure of Figures 4 and 5 is different from that of Figure 3 in respect of the total number of equations and parameters, since the number of time-dependent currents is not the same. The increase in the number of currents in the generic model of Figures 4 and 5 was necessitated by the requirement that the optimised model fit three AP datasets simultaneously, as opposed to the two AP datasets of Figure 3.

tab5
Table 5: Parameter values for random-paced LA optimisation.
tab6
Table 6: Initial variable values for random-paced optimisation.
fig5
Figure 5: Reconstructed membrane currents of the optimised model in Figure 4. From top to bottom: membrane currents in response to stimulation at pacing intervals (PIs) of 400, 200 ms and a random sequence of PIs.
3.2.3. Multidataset Optimisation with Combined Shared and Data-Specific Parameters

(i) Drug-Specific Case. The generic ionic model was simultaneously fitted to two spontaneous peripheral sinoatrial node (pSAN) AP datasets from Kodama et al. [21]: a control set of APs and the AP response following the application of E-4031, a specific blocker of rapid delayed-rectifier potassium channels. To fit the generic model to both datasets, four time-dependent currents and one leakage current were required. Current was chosen to correspond to . Furthermore, it was assumed that E-4031 acts to only alter the maximum conductance of channels, without modulating their kinetics. Therefore during optimisation, only one pharmacological-specific parameter , the maximum conductance of the , was optimised to have a distinct value for each of the two datasets. All other model parameter values were shared between the datasets, since any AP waveshape variation was assumed to be due to the blocking effect of E-4031 on alone. The optimised pSAN model was spontaneously active and AP fits to both datasets are shown in Figure 6, with RMS error between the model and corresponding data being 2.79 mV (control) and 3.91 mV (E-4031). Optimised parameter and initial values of model variables are given in Tables 7 and 8. Under control conditions, the optimised value was 2438.40 μS/cm2 but was reduced to 1562.28 μScm−2 under the action of E-4031. Model-generated membrane currents corresponding to these fits are shown in Figures 7(a) and 7(b). The amplitude of in the E-4031 model is evidently reduced from a mean peak value of 19.26 μA/cm2 to 12.73 μA/cm2, resulting in a prolongation of repolarization and a decreased maximum diastolic potential. Note that membrane currents , , , and , which are not directly affected by E-4031, still exhibit differences between the two simulations due to the voltage dependency of each current and its corresponding interaction with .

tab7
Table 7: Parameter values for drug-specific optimisation.
tab8
Table 8: Initial variable values for drug-specific optimisation.
fig6
Figure 6: Drug-specific multiobjective optimisation. Top and bottom panels illustrate the optimised model (solid traces) overlaid with experimental AP data (dashed traces) of peripheral sinoatrial node (pSAN) APs under control conditions (a) and in the presence of 0.1 μM E-4031, a selective blocker of channels (b).
fig7
Figure 7: Reconstructed membrane currents for the drug-specific optimisation. Four time-dependent currents and one leakage current were included in the model. , which was partially blocked by E-4031, is shown as a thick red dashed line.

(ii) Tissue-Specific Case. Spatial heterogeneity in AP waveforms is evident throughout the atria, likely due to the differential expression of ion channels from sinoatrial node (SAN) to atrial regions [19]. APs were recorded from central sinoatrial node (cSAN) and right atrial (RA) tissue-intact myocytes in one rabbit SAN preparation ( for each cell type). For this case of tissue-specific optimisation, parameters for the maximum conductance of each membrane current (an indicator of ion channel density) were set as data-specific parameters. Each cell type possessed distinct values for the , while all remaining parameters were shared: the assumption being that ion channels with the same kinetic properties are present throughout the whole tissue.

Figure 8 shows the fitted APs for both regions. The cSAN model was spontaneously active, whilst the atrial model was stimulated with rectangular pulses of 2 ms duration and 22 μA/cm2 amplitude. A total of seven time-dependent currents and one leakage current were required to fit both datasets simultaneously. The RMS errors between the optimised models and corresponding data were 1.78 mV (cSAN) and 2.48 mV (RA). Values of all optimised parameters and initial variables are given in Tables 9 and 10. Significant channel density differences (maximum conductances) exist between cSAN and RA, which contribute to the difference in AP waveform. Comparison of cSAN and RA membrane currents is plotted in Figure 9, with a separate panel for each current from both cell types. Interestingly, ionic currents with the same kinetic parameters can display a very different time course for each cell type, due to the direct effect of the differences in and the indirect effects of voltage dependency in each current.

tab9
Table 9: Parameter values for tissue-specific optimisation.
tab10
Table 10: Initial variable values for tissue-specific optimisation.
fig8
Figure 8: Tissue-specific multiobjective optimisation. Top and bottom panels illustrate the optimised model (solid traces) overlaid with experimental AP data (dashed traces), representing central sinoatrial node (cSAN) and right atrial (RA) APs recorded from the same tissue preparation.
fig9
Figure 9: Model-generated membrane currents corresponding to the tissue-specific optimisation of Figure 8. Seven time-dependent currents and one leakage current were included in the model. All model parameters were shared between the two cell types, with the exception of the maximum membrane conductance for each of the ionic currents (i.e., and , a total of eight parameters).

4. Discussion

In this study, a generic ionic model was optimised using a custom curvilinear gradient algorithm to fit a range of cardiac APs. The model was fitted to either single AP traces or simultaneously to multiple AP waveforms recorded from tissue-intact myocytes under different experimental conditions in rabbit SAN and/or atrial tissue preparations. AP waveforms could be well reproduced by the generic model, whose complexity lies intermediate between simple phenomenological formulations and biophysically detailed ionic models. Complex experimental data could be reproduced by the addition of extra ionic currents into the model. A major improvement over existing modelling approaches is that model parameters have been adjusted to accurately reproduce AP waveforms recorded under different pacing or pharmacological conditions from the same myocyte, reproducing complex AP characteristics while retaining physiologically realistic membrane current waveforms.

Membrane current kinetics of the generic model were expressed in terms of two first-order gates ( and ). It is possible to incorporate more complex activation kinetics such as raising the gating variables to powers greater than one. Such a modification would, in principle, allow fits to sigmoidal time courses of activation during voltage clamps: a property of many membrane currents [22]. Surprisingly, our fits to multiple AP data did not require sigmoidal kinetics for the membrane currents, as would be the case if we were to reproduce voltage-clamp data. However if desired, voltage-clamp data could be included as an additional dataset to be fitted alongside the other AP records simultaneously. Even if the data required membrane currents to reproduce a sigmoidal onset of activation, this could still be achieved with the simplified kinetic structure of our genetic model. For example, we have fitted the Hodgkin and Huxley (HH) [5] ( kinetics) and (m3 h kinetics) membrane currents in response to a voltage step from −60 mV to +40 mV, using a total of two generic currents (for ) and four generic currents (for ) respectively (not shown), indicating that even with highly simplified kinetics, the generic model is still able to reproduce a wide range of experimental data. The modeller can decide whether to amalgamate any generic currents obtained with our method into more complex formulations, as a first step towards formulating a more biophysically detailed model if desired.

Compared with previous studies using simplified cell formulations such as the Fitzhugh-Nagumo or Fenton-Cherry models [4, 14, 23], our generic model could accurately reproduce spontaneous AP waveforms recorded from cSAN, pSAN as well as paced AP waveforms from RA and LA myocytes. Although previously published phenomenological models were able to reproduce restitution curves and reasonable AP waveforms, we regard it necessary to reproduce accurate AP morphologies in modelling electrophysiological dynamics [3, 24, 25]. With our approach, a user-defined number of membrane currents can be defined, providing higher flexibility in reproducing even more complex electrophysiological dynamics. In order to retain the simplified nature of the model, additional ionic currents were included only when necessary (i.e., the target RMS error could not be achieved). Moreover, because of the similar formulation of each ionic current, many of these can be recombined if they are found to follow similar time course profiles during optimisation, facilitating the process of model reduction. It is important to note that prior to optimisation, we make no assumptions as to the ionic identity of each current, but allow the fitting process to determine the current density and kinetics based solely on the AP data. The only exception was for the drug-specific scenario (Figures 6 and 7), where, prior to optimisation, the maximum conductance of current was chosen as the only parameter to be altered by E-4031, effectively preidentifying this current with for this optimisation run only. We could equally have chosen any other conductance parameter in the model. Although is preidentified with in this case, does not necessarily correspond to this membrane current for the other optimisations results of this study. The generic model approach outlined here provides a promising tool for tissue or whole-heart simulations due to its simplified nature and, therefore, computational efficiency.

Our multidataset-based optimisation results suggest that multiobjective experimental AP data can improve ionic model optimisation and predictive power of the model, provided the additional data includes information not present in the original dataset(s). Model optimisation using multiple data with similar AP characteristics, such as AP duration or amplitude, will tend to simply maintain the number and location of local minima on the objective surface. However, introducing additional datasets with extra information will smooth the objective surface by reducing the amplitude of any “surface ripple,” since each parameter point on the surface must now fit multiple data (see Figure 10). In addition, new datasets with distinct “information” will introduce more local minima onto the objective surface, making its topology more complex and thus confounding the search for the global optimum. We found that for multiple data optimisation, it was much more challenging to fit all datasets simultaneously, particularly if there were stringent constraints on parameters and parameter values were shared between the datasets. At the same time, the credibility of the ionic model is enhanced by its ability to simultaneously reproduce data obtained under variable experimental conditions

706195.fig.0010
Figure 10: Normalised objective surface reproduced by the RMS error between model and experimental APs, against two optimised parameters and representing reversal potential and maximum membrane conductance . Upper plot, 2D objective surface of multi-dataset based optimisation. Note the “noisy” local area indicated by the red arrow. The global minimum is labelled by the black arrow. Lower plot, 2D objective surface of single-dataset based optimisation. Inset: Zoomed 1D local objective surface near global minimum. When optimising the model to fit multiple datasets simultaneously, the number of local minima will be increased. However, each local minimum will be shallower compared to those of the single-dataset based surface.

In addition, reconstructed membrane currents of the optimised models were found to follow physiologically realistic waveforms in SAN and atrial rabbit myocytes, consistent with current profiles predicted in existing biophysically detailed models [26, 27]. It is important to note that membrane currents with unrealistic time course profiles may perfectly reproduce a single AP record dataset, but these will generally fail to simultaneously reproduce multiple experimental data. Our results suggest that the multiobjective fitting approach can be used to accurately reconstruct underlying membrane current dynamics, since the additional information provided by the multiple data was important for their accurate reconstruction.

From the experimental APs recorded in response to pacing at uniform and random intervals, it can be seen that for both the PI = 200 ms and random paced datasets, there was some beat-beat variation in the time course of underlying membrane currents (Figure 5). At high pacing frequencies, a second AP is elicited shortly after the previous one, having a reduced AP duration and refractory period due to the reduction in magnitude of inward ionic currents. The generic model was optimised to simultaneously fit AP responses to regular and random pacing and could therefore reproduce this AP alternans at higher rates, which is important in understanding and simulating the pathophysiological mechanisms underlying reentrant activity and electrical remodelling in atrial fibrillation [28].

Many electrophysiological studies have employed dynamic AP clamp in the presence of selective ion channel blockers to visualise the time course of ionic currents underlying the AP. However, these methods suffer from the limitation that it is not currently possible to simultaneously visualise all membrane currents present in a given myocyte. An alternative approach may be to use integrative ionic model simulations which can reproduce cellular electrophysiological behaviour under multiple experimental conditions. The generic model of this study was optimised to fit multiple AP records recorded during control conditions and also under the influence of a selective ion channel blocker (E-4031). It was assumed that in the two conditions, all model parameters shared identical values, with the exception of the maximum membrane conductance of the blocked current. We have found that fitting to a single AP dataset does not provide unique reconstructions of the time course of underlying membrane currents. Nevertheless, the membrane currents reconstructed using the multiobjective drug-specific optimisation of the present study reveal remarkably similar time course profiles to other existing biophysically detailed models. For example, rabbit SAN [29, 30] and rabbit atrial [26, 30] single cell ionic models display a similar time course to our , suggesting that this optimised current has been appropriately constrained by the drug-specific data optimisation. In contrast, some experimental AP-clamp data [31, 32] indicate a slower onset of , suggesting that we have not accurately isolated the time course of this current component, perhaps through inclusion of an component. Additional datasets for use with multiobjective optimisation, such as, for example, voltage-clamp data in the presence of the drug, may further help in more accurately identifying the kinetics of this current.

The generic model approach outlined here could also be used to simulate the effects of antiarrhythmic drugs on whole heart or tissue simulations, due to the computational simplicity of the model. The approach could even be used with pharmacological agents known to have multiple effects, such as the partial block of more than one current. The user would only need to specify which model parameters are shared between datasets and which are drug-specific.

The tissue-specific optimisation approach also allowed the generic model to be fitted to heterogeneous APs from different myocyte types in the same tissue preparation. It was assumed in this case that all model parameters were shared between datasets, except the maximum membrane conductance of each current. As was the case for the drug-specific optimisation, the reconstructed generic currents in the model resembled known profiles of ionic currents, based on their similarity to current profiles obtained from biophysically detailed models of cSAN and atrial rabbit myocytes, as well as published data of pharmacologically isolated ionic currents obtained using the dynamic AP clamp technique [3338] (see Figure 11). In addition, comparison of the magnitudes of and ( and , resp.) between SAN and atrial myocytes is in agreement with experimental data from isolated right atrial preparations [39, 40].

fig11
Figure 11: Comparison of the time course of generic model currents with experimental AP clamp data and biophysically-detailed model simulations.

Furthermore, the excellent AP fits obtained for the tissue-specific case (Figure 8) indicate that this multiobjective optimisation approach could also provide an extension to previously published gradient models of cardiac tissue electrical activity, particularly models of SAN-atrial interaction. Compared with Lovell et al.’s work [29], all kinetic model parameters were shared by each dataset in this study; therefore any spatial regional differences are only reflected through variation in maximal membrane conductances. Compared to the gradient model published by Zhang et al. [30], the tissue-specific model of this study could accurately fit SAN and RA APs by optimising shared kinetic parameters. Compared with the work of Syed et al. [19] and Dastgheib et al. [18], who only optimised maximum ion channel conductance parameters [18, 19], all model parameters were included in the tissue-specific optimisation of Figure 8. We believe that estimating only conductance parameters while fixing ion channel kinetics parameters will excessively limit the parameter search space, reducing the accuracy of model fits, particularly when optimising against multiple datasets.

Like most existing models, our generic model has certain limitations, representing compromises which are necessary to achieve simplified and computationally efficient descriptions of membrane current kinetics. While the model can accurately fit multiple AP data, we have not incorporated intracellular calcium cycling, changes in intracellular ion concentrations, metabolites, or ionic pumps and exchangers. These changes, if present, would impact the AP waveshape through our reconstructed membrane currents, which we assume to simply consist of two first-order voltage-dependent ( and ) gating processes. It is also important to note that optimisation of nonlinear models, particularly those with large numbers of parameters, may lead to nonunique parameter estimates. The generic model of this study is no exception. The symmetrical nature of the membrane current formulations indicates that the parameters of any two membrane currents can be interchanged to produce identical AP waveforms, due to the identical formulations for each current. A similar argument would hold for the gating variables: their kinetic parameters can be interchanged within any current due to their symmetrical formulations. From such simple considerations, it can be concluded that parameters of the generic model cannot be uniquely determined, unless nonsymmetrical upper/lower bounds are imposed on individual parameters. However, we found the model currents can converge to similar waveforms, regardless of the initial parameter values used. These results indicate that the additional information provided by the multiple data was important for accurate reconstruction of membrane currents. In general, this was not possible when fitting the model to single datasets, even though the AP record itself could be well fitted [20, 41]. In other words, model behaviours, as opposed to parameters, tend to be well identified when extra experimental information is incorporated into the optimisation.

5. Conclusion

We have presented a generic model of cardiac electrical activity, capable of accurately reproducing action potential waveforms from multiple experimental data in any given myocyte. Furthermore, we have shown that our generic ionic model and multiobjective optimisation approach described in this study can provide an effective and efficient means to reconstruct the profiles of hidden ionic currents underlying the AP. Multiobjective fitting to multiple AP datasets appears to provide stringent constraints on the dynamics of underlying membrane currents, yielding reconstructed membrane current time course profiles in agreement with existing studies. The generic approach will allow the efficient computation of complex electrophysiological dynamics in whole heart simulations and will provide a valuable tool in elucidating the ionic mechanisms underlying cardiac electrical activity.

References

  1. D. C. Sigg, P. A. Iaizzo, Y. F. Xiao, and B. He, Cardiac Electrophysiology Methods and Models, Springer, New York, NY, USA, 2010.
  2. R. L. Winslow, D. F. Scollan, A. Holmes, C. K. Yung, J. Zhang, and M. S. Jafri, “Electrophysiological modeling of cardiac ventricular function: from cell to organ,” Annual Review of Biomedical Engineering, vol. 2, no. 2000, pp. 119–155, 2000. View at Scopus
  3. M. Fink, S. A. Niederer, E. M. Cherry et al., “Cardiac cell modelling: observations from the heart of the cardiac physiome project,” Progress in Biophysics and Molecular Biology, vol. 104, no. 1–3, pp. 2–21, 2011. View at Publisher · View at Google Scholar · View at Scopus
  4. R. Fitzhugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical Journal, vol. 1, no. 6, pp. 445–466, 1961. View at Publisher · View at Google Scholar
  5. A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” The Journal of Physiology, vol. 117, no. 4, pp. 500–544, 1952. View at Scopus
  6. M. Courtemanche, W. Skaggs, and A. T. Winfree, “Stable three-dimensional action potential circulation in the Fitzhugh-Nagumo model,” Physica D, vol. 41, no. 2, pp. 173–182, 1990. View at Scopus
  7. B. Y. Kogan, W. J. Karplus, B. S. Billett, A. T. Pang, H. S. Karagueuzian, and S. S. Khan, “The simplified FitzHugh-Nagumo model with action potential duration restitution: effects on 2D wave propagation,” Physica D, vol. 50, no. 3, pp. 327–340, 1991. View at Scopus
  8. R. R. Aliev and A. V. Panfilov, “A simple two-variable model of cardiac excitation,” Chaos, Solitons and Fractals, vol. 7, no. 3, pp. 293–301, 1996. View at Publisher · View at Google Scholar · View at Scopus
  9. M. P. Nash and A. V. Panfilov, “Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias,” Progress in Biophysics and Molecular Biology, vol. 85, no. 2-3, pp. 501–522, 2004. View at Publisher · View at Google Scholar · View at Scopus
  10. C. H. Luo and Y. Rudy, “A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction,” Circulation Research, vol. 68, no. 6, pp. 1501–1526, 1991. View at Scopus
  11. G. W. Beeler and H. Reuter, “Reconstruction of the action potential of ventricular myocardial fibres,” The Journal of Physiology, vol. 268, no. 1, pp. 177–210, 1977. View at Scopus
  12. F. H. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation,” Chaos, vol. 8, no. 1, pp. 20–47, 1998. View at Publisher · View at Google Scholar
  13. F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans, “Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity,” Chaos, vol. 12, no. 3, pp. 852–892, 2002. View at Publisher · View at Google Scholar · View at Scopus
  14. E. M. Cherry, J. R. Ehrlich, S. Nattel, and F. H. Fenton, “Pulmonary vein reentry—properties and size matter: insights from a computational analysis,” Heart Rhythm, vol. 4, no. 12, pp. 1553–1562, 2007. View at Publisher · View at Google Scholar · View at Scopus
  15. S. Nattel, B. Burstein, and D. Dobrev, “Atrial remodeling and atrial fibrillation: mechanisms and implications,” Circulation, vol. 1, no. 1, pp. 62–73, 2008. View at Publisher · View at Google Scholar · View at Scopus
  16. M. Courtemanche, R. J. Ramirez, and S. Nattel, “Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model,” The American Journal of Physiology, vol. 275, no. 1, pp. H301–H321, 1998. View at Scopus
  17. E. M. Cherry and F. H. Fenton, “A tale of two dogs: analyzing two models of canine ventricular electrophysiology,” The American Journal of Physiology, vol. 292, no. 1, pp. H43–H55, 2007. View at Publisher · View at Google Scholar · View at Scopus
  18. Z. S. Dastgheib, A. Azemi, M. Khademi et al., “Identification of ionic conductances in a reentry model of ventricular myocardium cells,” Journal of Applied Sciences, vol. 9, no. 3, pp. 555–560, 2009. View at Scopus
  19. Z. Syed, E. Vigmond, S. Nattel, and L. J. Leon, “Atrial cell action potential parameter fitting using genetic algorithms,” Medical and Biological Engineering and Computing, vol. 43, no. 5, pp. 561–571, 2005. View at Publisher · View at Google Scholar · View at Scopus
  20. S. Dokos and N. H. Lovell, “Parameter estimation in cardiac ionic models,” Progress in Biophysics and Molecular Biology, vol. 85, no. 2-3, pp. 407–431, 2004. View at Publisher · View at Google Scholar · View at Scopus
  21. I. Kodama, M. R. Boyett, M. R. Nikmaram, M. Yamamoto, H. Honjo, and R. Niwa, “Regional differences in effects of E-4031 within the sinoatrial node,” The American Journal of Physiology, vol. 276, no. 3, pp. H793–H802, 1999. View at Scopus
  22. B. Hille, Ion Channels of Excitable Membranes, Sinauer Associates, Sunderland, Mass, USA, 3rd edition, 2001.
  23. A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton, “Minimal model for human ventricular action potentials in tissue,” Journal of Theoretical Biology, vol. 253, no. 3, pp. 544–560, 2008. View at Publisher · View at Google Scholar · View at Scopus
  24. E. M. Cherry and F. H. Fenton, “Suppression of alternans and conduction blocks despite steep APD restitution: electrotonic, memory, and conduction velocity restitution effects,” The American Journal of Physiology, vol. 286, no. 6, pp. H2332–H2341, 2004. View at Publisher · View at Google Scholar · View at Scopus
  25. K. H. W. J. T. Tusscher and A. V. Panfilov, “Alternans and spiral breakup in a human ventricular tissue model,” The American Journal of Physiology, vol. 291, no. 3, pp. H1088–H1100, 2006. View at Publisher · View at Google Scholar · View at Scopus
  26. D. S. Lindblad, C. R. Murphey, J. W. Clark, and W. R. Giles, “A model of the action potential and underlying membrane currents in a rabbit atrial cell,” The American Journal of Physiology, vol. 271, no. 4, pp. H1666–H1696, 1996. View at Scopus
  27. S. Dokos, B. Celler, and N. Lovell, “Ion currents underlying sinoatrial node pacemaker activity: a new single cell mathematical model,” Journal of Theoretical Biology, vol. 181, no. 3, pp. 245–272, 1996. View at Publisher · View at Google Scholar · View at Scopus
  28. S. Nattel, “New ideas about atrial fibrillation 50 years on,” Nature, vol. 415, no. 6868, pp. 219–226, 2002. View at Publisher · View at Google Scholar · View at Scopus
  29. N. H. Lovell, S. L. Cloherty, B. G. Celler, and S. Dokos, “A gradient model of cardiac pacemaker myocytes,” Progress in Biophysics and Molecular Biology, vol. 85, no. 2-3, pp. 301–323, 2004. View at Publisher · View at Google Scholar · View at Scopus
  30. H. Zhang, A. V. Holden, I. Kodama et al., “Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node,” The American Journal of Physiology, vol. 279, no. 1, pp. H397–H421, 2000.
  31. G. M. Faber, J. Silva, L. Livshitz, and Y. Rudy, “Kinetic properties of the cardiac L-type Ca2+ channel and its role in myocyte electrophysiology: a theoretical investigation,” Biophysical Journal, vol. 92, no. 5, pp. 1522–1543, 2007. View at Publisher · View at Google Scholar · View at Scopus
  32. J. C. Hancox, A. J. Levi, and H. J. Witchel, “Time course and voltage dependence of expressed HERG current compared with native “rapid” delayed rectifier K current during the cardiac ventricular action potential,” Pflugers Archiv—European Journal of Physiology, vol. 436, no. 6, pp. 843–853, 1998. View at Publisher · View at Google Scholar · View at Scopus
  33. T. Doerr, R. Denger, and W. Trautwein, “Calcium currents in single SA nodal cells of the rabbit heart studied with action potential clamp,” Pflugers Archiv—European Journal of Physiology, vol. 413, no. 6, pp. 599–603, 1989. View at Scopus
  34. N. Kim, M. B. Cannell, and P. J. Hunter, “Changes in the calcium current among different transmural regions contributes to action potential heterogeneity in rat heart,” Progress in Biophysics and Molecular Biology, vol. 103, no. 1, pp. 28–34, 2010. View at Publisher · View at Google Scholar · View at Scopus
  35. L. Yue, J. Feng, R. Gaspo, G. R. Li, Z. Wang, and S. Nattel, “Ionic remodeling underlying action potential changes in a canine model of atrial fibrillation,” Circulation Research, vol. 81, no. 4, pp. 512–525, 1997. View at Scopus
  36. M. E. Mangoni, B. Couette, L. Marger, E. Bourinet, J. Striessnig, and J. Nargeot, “Voltage-dependent calcium channels and cardiac pacemaker activity: from ionic currents to genes,” Progress in Biophysics and Molecular Biology, vol. 90, no. 1–3, pp. 38–63, 2006. View at Publisher · View at Google Scholar · View at Scopus
  37. G. Berecki, J. G. Zegers, R. Wilders, and A. C. van Ginneken, “Cardiac channelopathies studied with the dynamic action potential-clamp technique,” Methods in Molecular Biology, vol. 403, pp. 233–250, 2007. View at Scopus
  38. A. Nygren, C. Fiset, L. Firek et al., “Mathematical model of an adult human atrial cell: the role of K+ currents in repolarization,” Circulation Research, vol. 82, no. 1, pp. 63–81, 1998. View at Scopus
  39. H. Honjo, M. R. Boyett, I. Kodama, and J. Toyama, “Correlation between electrical activity and the size of rabbit sino-atrial node cells,” The Journal of Physiology, vol. 496, no. 3, pp. 795–808, 1996. View at Scopus
  40. M. Lei and M. R. Boyett, “Inhibition of transient outward current, it(to), by flecainide and quinidine in rabbit isolated sinoatrial node cells,” The Journal of Physiology, vol. 511, pp. 78–79, 1998.
  41. A. X. Sarkar and E. A. Sobie, “Regression analysis for constraining free parameters in electrophysiological models of cardiac cells,” PLoS Computational Biology, vol. 6, no. 9, Article ID e1000914, 2010. View at Publisher · View at Google Scholar · View at Scopus