Research Article  Open Access
Optimisation of a Generic Ionic Model of Cardiac Myocyte Electrical Activity
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 userdefined number of ionic currents, employing twogate HodgkinHuxley 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 voltageclamp protocols, AP membrane current dynamics, and Ca^{2+} alternans [2, 3]. However, the everincreasing number of variables in such models renders them computationally expensive when integrated into higherdimensional 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 FitzhughNagumo (FHN) formulation published in 1961 [4], simplifying the fourvariable HodgkinHuxley (HH) nerve axon model [5] into a twovariable formulation by eliminating gating variables with rapid time constants. Since its publication, FHNtype models have been frequently used in multicellular tissue simulations [6–9]. 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 threevariable FentonKarma 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 FentonCherry 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 everincreasing complexity of such models makes parameter estimation a highly difficult and timeconsuming 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 largescale optimisation of cardiac cell ionic models.
In response to these challenges, we have formulated a simplified HHtype generic model to accurately reconstruct a range of AP waveforms recorded from tissueintact 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 highfrequency paced stimulation. Model parameters were estimated using a computationally simple, multiobjective AP optimisation approach based on a custom curvilineargradient 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/cm^{2}) is the specific membrane capacitance, (nA/cm^{2}) refers to the timeindependent (leakage) current, (nA/cm^{2}) is the applied stimulus, (nA/cm^{2}) denotes the th timedependent ionic current density, and is the userdefined number of timedependent 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 timedependent currents are represented by a twogate HH scheme: where (μS/cm^{2}) is the maximum ion channel conductance, and are dimensionless gating variables, and (mV) is the reversal potential for membrane current . The timeindependent leakage current is described by where (μS/cm^{2}) 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 gradientbased leastsquares optimisation method, combining the advantages of both Newton and steepestdescent 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 leastsquares 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. Dataspecific 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 dataspecific parameters. This is the case for the drugspecific and tissuespecific optimisation scenarios presented. For dataspecific 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. SANRA or LA appendage tissue preparations were dissected and placed in a recording chamber, superfused with Tyrode’s solution and oxygenated with 95% O_{2} and 5% CO_{2} 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 USB6251 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 HHtype ionic model with two timedependent 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 tissueintact 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 modelgenerated 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/cm^{2} 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.


(a)
(b)
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 threecurrent 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 peaktopeak 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.
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 timedependent 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/cm^{2} in amplitude were used to elicit APs in the model. For a PI of 200 ms, the model generated AP characteristics and corresponding timedependent ionic currents revealed beattobeat 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).
(a)
(b)
To test the reliability of the previous multipledataset 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 modelgenerated 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 timedependent ionic current waveforms are nearly identical and physiologically reasonable, except for differences in their scaling (Figure 3(b)).


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 beatbeat variations in AP morphology were demonstrated in the randomly paced dataset, as shown in Figure 4 (lower panel).
(a)
(b)
(c)
Using the generic ionic model, a total of seven timedependent 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/cm^{2} 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 randompaced 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 timedependent 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.


(a)
(b)
(c)
3.2.3. Multidataset Optimisation with Combined Shared and DataSpecific Parameters
(i) DrugSpecific 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 E4031, a specific blocker of rapid delayedrectifier potassium channels. To fit the generic model to both datasets, four timedependent currents and one leakage current were required. Current was chosen to correspond to . Furthermore, it was assumed that E4031 acts to only alter the maximum conductance of channels, without modulating their kinetics. Therefore during optimisation, only one pharmacologicalspecific 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 E4031 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 (E4031). 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/cm^{2} but was reduced to 1562.28 μScm^{−2} under the action of E4031. Modelgenerated membrane currents corresponding to these fits are shown in Figures 7(a) and 7(b). The amplitude of in the E4031 model is evidently reduced from a mean peak value of 19.26 μA/cm^{2} to 12.73 μA/cm^{2}, resulting in a prolongation of repolarization and a decreased maximum diastolic potential. Note that membrane currents , , , and , which are not directly affected by E4031, still exhibit differences between the two simulations due to the voltage dependency of each current and its corresponding interaction with .
