The traditional cardiac model-building paradigm involves constructing a composite model using data collected from many cells. Equations are derived for each relevant cellular component (e. g. , ion channel, exchanger) independently. After the equations for all components are combined to form the composite model, a subset of parameters is tuned, often arbitrarily and by hand, until the model output matches a target objective, such as an action potential. Unfortunately, such models often fail to accurately simulate behavior that is dynamically dissimilar (e. g. , arrhythmia) to the simple target objective to which the model was fit. In this study, we develop a new approach in which data are collected via a series of complex electrophysiology protocols from single cardiac myocytes and then used to tune model parameters via a parallel fitting method known as a genetic algorithm (GA). The dynamical complexity of the electrophysiological data, which can only be fit by an automated method such as a GA, leads to more accurately parameterized models that can simulate rich cardiac dynamics. The feasibility of the method is first validated computationally, after which it is used to develop models of isolated guinea pig ventricular myocytes that simulate the electrophysiological dynamics significantly better than does a standard guinea pig model. In addition to improving model fidelity generally, this approach can be used to generate a cell-specific model. By so doing, the approach may be useful in applications ranging from studying the implications of cell-to-cell variability to the prediction of intersubject differences in response to pharmacological treatment.
Mathematical models of cardiac electrophysiology trace their roots to Hodgkin and Huxley’s seminal work from 1952 [1]. Since then many models have been developed describing cardiac electrophysiology in a number of species and cell types helping to make modeling an integral part of cardiac research [2–5]. The typical method for model development and parameterization is a bottom-up approach. Individual ionic membrane currents are characterized using voltage-clamp experiments from which mathematical equations are derived [6,7]. Although it has led to many advances, this traditional approach to model development has several limitations, including: Several studies have looked into how to improve model parameterization. Approaches in cardiac myocyte modeling have included the parameterization of individual channel dynamics, typically when making more complex Markov models of ionic currents [6,15–17]. Whole-cell optimization approaches have focused on generating models that can match action potentials from different types of cardiomyocytes, using both simple models consisting of a few generic currents [13,18–20] and more physiologically detailed ionic models [21–25]. A resulting synthesis is that optimization results are improved when models are fit to data beyond a single action potential, e. g. , action potentials from multiple pacing rates [13,19,22] or voltage waveforms during varying current injection [20,21,24]. In particular, using a global search heuristic applied to an ionic model, Syed et al. demonstrated that it is feasible to estimate conductance parameters for experimental data and showed that the fits improved when using data recorded during multiple periodic pacing frequencies [22]. Sarkar and Sobie presented a much simpler, but more approximate, linear regression strategy to estimate model conductances based on biomarkers from simulated model output and have used it to investigate how specific conductances relate to particular model outputs [26]. In neuroscience, considerably more research has been carried out on parameter estimation problems (e. g. , [27–30]) and a few studies have developed protocols that allow parameterization of cell-specific models [31,32]. However, these protocols are not directly applicable to cardiac myocytes, due to intrinsic differences in electrophysiological behavior between neurons and cardiomyocytes. Here, we present a novel strategy to develop cardiac models by pairing dynamically rich electrophysiology protocols with powerful computational parameter fitting methods. We first developed novel electrophysiology protocols that probe the dynamics of a subset of ionic currents in an intact cardiac myocyte without ion channel inhibitors, agonists, or unphysiological ion concentrations. The protocol consists of (1) stochastic current-clamp stimulation and (2) multi-step voltage clamping. As will be discussed, stochastic stimulation represents a quick method to sample the rate-dependent cardiac dynamics, while the multi-step voltage-clamp protocol is designed to highlight individual currents using a tailored sequence of holding potentials. Based on the assumption that ion-channel kinetics are preserved among (healthy) subjects while maximal conductances vary as a result of differences in expression levels, the resulting data are used to estimate maximal conductance values of several ionic currents and maximal flux of calcium ion transporters in the model. Because of the complexity of the data, hand parameter tuning is not feasible; thus, we utilized a genetic algorithm (GA), which is an efficient method for such complex optimization applications [33]. The approach was first developed and validated computationally. It was then used to develop cell-specific models of isolated guinea pig ventricular myocytes that simulate the electrophysiological dynamics significantly better than does a standard guinea pig model.
We first developed our parameter estimation strategy using a guinea pig ventricular myocyte model (Faber and Rudy [34], the “FR” model) and tested the ability of the optimization procedure to return the original parameter values. Traditionally, one of the main criteria for cardiac electrophysiology model quality is the ability of a model to describe the cardiac action potential. Therefore, we first ran the parameter estimation using a single FR model action potential as the target objective. Nine model parameters, describing maximal conductances of ionic currents [the sodium current (INa), the L-type calcium current (ICaL), the T-type calcium current (ICaT), the inwardly rectifying potassium (IK1), the rapid and slow delayed rectifier potassium currents (IKr and IKs), the plateau potassium current (IKp), and the sarcolemmal calcium pump current (IpCa) ] and the maximal flux of the sarcoplasmic reticulum Ca2+-ATPase (JSERCA) were estimated using a GA technique. A GA optimization is initialized by a population of models with different parameter values. We used a population size of 500 model instantiations generated by randomly drawing values for the nine parameters from a range of 0. 01–299% of the published value. The GA methodology uses ideas from evolution [23,33]. In GA terminology, the initial state is referred to as generation 0 and the 500 models as individuals. Fig 1 A–1 C shows three different individuals in generation 0. Most of these generate action potentials that are very different from the target action potential (Fig 1A and 1C). However, by chance, a few of the 500 individuals provide reasonably good fits to the action potential, even though the parameter values are very different from those of the baseline model (parameter scaling of 1, Fig 1B). The optimization proceeds in generations (steps), for which the GA applies crossover (parameter swapping), mutation (parameter variation), and selection (discarding poorly performing models) to increase the fitness of the population (reduce the error between model output and target objective). We ran the GA for 100 generations, as this was sufficient for the error of both the population average and the best individual to reach a minimum (Fig 1 F and 1G). Although the optimized model action potential matches the optimization objective to a very high degree, the estimated parameter set does not match that of the FR model (Fig 1E). This is consistent with previous results showing that if only a single action potential is used for parameterization, cardiac models may be overdetermined and more than a single set of parameter values can describe that action potential [20,21,26]. The duration of the action potential, and to a lesser extent its morphology, varies with stimulation interval and history. Models tuned to single action potentials during periodic pacing (as in Fig 1) would not be expected to accurately reproduce such dynamics. A more complete method to probe cellular behavior is the restitution portrait [35], which is a systematic, but prohibitively long, mapping of this rate dependence. An alternative approach that would significantly reduce the protocol duration, while maintaining some dynamic information, is random sampling of the rate dependence. To accomplish this, we utilized a stochastic stimulation protocol containing 11 randomly timed stimuli delivered over 5 s. When applied to the FR model, the stochastic stimulation results in considerable action potential variability (Fig 2 A and 2B). We used this stochastic stimulation protocol and resulting voltage response as an optimization sequence to test the extent to which dynamic stimulus timing would improve the parameter estimation. Because the GA parameter estimation is a stochastic method, it was run 10 times with 10 different initial populations. For each run, we selected the best individual, i. e. , the model instantiation with the best fit to the target voltage trace. All 10 best individuals matched the voltage trace very closely (Fig 2 A and 2B show the best one) as was the case for the single action potential fitting. Compared to the single action potential, the stochastic stimulation leads to a modest overall improvement of the parameter estimation, but it did notably better in determining the maximal conductances of IKr, ICaL, and IKs (Fig 2C). However, as shown in Fig 2C, a few parameters remain incorrectly estimated (maximal conductances of ICaL and IKs), and some are estimated with a large spread (maximal conductance values of IKp, IpCa, ICaT, and IKr). Sensitivity and correlation analyses illuminate why these parameters are less well determined—they have low sensitivity (in which case they have minimal effect on the fitting objective, making them difficult to probe) and/or they have inter-correlations (in which case two or more parameters lack independent contributions to the fitting objective and are therefore difficult to discriminate) (S1–S4 Figs and S1 Text). To measure the predictive capabilities of the optimized models, we presented the 10 best individuals from the GA runs with a novel stochastic stimulation sequence (S6 Fig). When subjected to this new prediction sequence, both the individuals trained with the stochastic stimulation sequence and those trained to the single action potential matched the baseline FR model response well (S6 Fig), but the stochastic stimulation protocol lead to the better match, as shown by the smaller prediction error (error between optimized model and target during the prediction sequence) in Fig 2D. Thus, although the parameter recovery seem only modestly improved for the stochastic stimulation compared to the single action potential target, the stochastic stimulation outperforms the single action potential in that it results in models that are significantly better at predicting the response to a novel set of stimuli. To improve parameter estimation accuracy, more improvement is typically gained from adding measurements of a different state variable than adding additional measurements of the same state variable [36]. This suggests that a longer stochastic stimulation protocol is unlikely to yield much improvement. This idea is in line with our finding above (Figs 2D and S6) that models optimized to the 5s stochastic stimulation protocol matched well when subjected to a novel stochastic stimulation sequence. Therefore, to improve the parameter estimation accuracy, we added a multi-step voltage clamp protocol to the objective function. Traditionally, voltage clamp is applied to a cell as a set of holding potentials varied systematically either in its timing or its potential to characterize one particular current. We developed a voltage clamp protocol consisting of a sequence of holding potentials, with each step designed to emphasize specific currents relative to the others (Fig 3). The rationale is that if a particular conductance contributes most of the total current for a particular holding potential, then only models fit with a correct value of that conductance will reproduce the current target for that phase of the protocol and therefore have a low fitting error. Our 6s long voltage clamp protocol effectively isolates IK1, ICaL, and IKs as shown by the disproportionally large contributions of these currents in step -120 mV, +20 mV, and +40 and -30 mV, respectively (Fig 3 B–3F). Therefore, we hypothesize that this protocol will directly improve the conductance estimation for these currents. Indeed, using the voltage clamp protocol as the objective during an optimization recovers the conductances for IK1 and IKs very accurately (Fig 4B). The estimation of the ICaL conductance is very close to 1, but is slightly overestimated in all runs. Less predictively, IpCa and IKp were also estimated more precisely than during stochastic pacing alone (Fig 4B). However, a few conductances were estimated poorly (in particular JSERCA and ICaT) and models optimized based on voltage clamp data alone were, not surprisingly, inferior at predicting complex action potential dynamics during stochastic pacing (Fig 4C). The extension of the target objective to include the multi-step voltage clamp protocol results in a joined unitless error function (Eq 3 in the Methods). As both the stochastic pacing recording and the voltage clamp data is fit increasingly well during optimization, the error contribution from each sequence decreases (Fig 4A, left). Although the main contribution to the total error comes from the stochastic pacing segment, the error from voltage clamp segment drops more during the optimization process, suggesting that both protocols help constrain the parameters. Running the optimization with the combined objective does indeed lead to improved accuracy of the parameter estimation, with all nine current parameters being recovered to within one standard deviation (orange symbols, Fig 4B). For six of the nine current parameters, the combined protocol results in parameters whose mean estimates are closer to 1 and/or have less variational spread than either of the individual protocols alone (Fig 4B). Only currents that were estimated very accurately by the voltage clamp protocol alone (IKs, IK1, and IpCa) did not show improvement with the combined protocol. For some of the currents, one protocol segment is clearly better than the other in terms of parameter recovery (e. g. , stochastic pacing for INa and voltage clamp for IKs). However, for other currents, in which both individual protocol segments result in off-target outcomes, the combined protocol produces estimates spanning 1 (e. g. , IKp). Such improvement is consistent with the combined protocol restraining parameter space and avoiding local minima. Again, we tested the ability of the 10 best individuals to predict the response to a stochastic stimulation sequence to which they were not fit. The prediction error of the 10 individuals from the combined objective function runs was significantly lower than the error for the individuals that were estimated using only the stochastic stimulation protocol (Fig 4C). Hence, the combined stochastic pacing and voltage clamp protocol improves both parameter recovery and prediction performance. To further improve the estimation results, the results of the first 10 GA runs were used as the new parameter bounds for a second set of runs (see Methods), e. g. 0. 01% to 299% changed to 92. 8–114. 9% for INa. Note that this method only works when the fits of the first 10 runs span the correct solution as is the case with the combined protocol. During this second, local, iteration, better fits are generated causing the error for both the voltage clamp and the current segments, as well as the total error for the best individual, to drop (Fig 4A, right). Thus, using this iterative approach, the error bounds around the estimated parameter values decreased (magenta symbols, Fig 4B) and the prediction error reduced markedly (Fig 4C). In summary, the combined protocol, consisting of stochastic stimulation and multi-step voltage clamp, allows accurate parallel estimation of eight maximal conductance values and maximal pump rate of SERCA for the FR guinea pig ventricular model. Such validation simulations laid the groundwork for using the method to fit computational models to real cardiac cell data. The parameter estimation method was next applied to four guinea pig left basal ventricular myocytes from four different animals. Each cell was subjected to the stochastic stimulation protocol in current clamp mode, followed by the multi-step voltage clamp protocol, using the perforated patch clamp technique. All four myocytes exhibited action potentials and membrane current responses that were very different from the baseline FR model (Fig 5 shows output from one cell, S7–S9 Figs presents the results from the remaining three cells). In particular, their action potentials were substantially longer than those of the FR model and their current response to prolonged depolarization was substantially smaller. For each cell, the GA estimate from the experimental data fit much better than did the FR model (Fig 5 and S7–S9 Figs). In particular, the optimization leads to very accurate voltage dynamics, which is important for arrhythmogenesis prediction. The total current is fit less well, potentially due to mismatch in ion channel kinetics (see Discussion). Overall, the optimization results in more accurate predictions, with the prediction error being an order of magnitude lower for the fitted models than for the FR model (Fig 5C and 5D). The dissimilarities between the original FR model and the experimental data led to considerable changes in the estimated values for the model parameters for all four cells (Fig 6). Interestingly, these changes were qualitatively similar between all four myocytes for most of the parameters, indicating conserved differences between our experimental data and the FR model. In particular, IKs and IKp are scaled down significantly and JSERCA is slightly reduced. In contrast, for all four cells, maximal conductance of IKr and IK1 are increased around 2-fold compared to the FR model, while ICaL is slightly increased. The results for INa show variation among cells, with a significant increase for three out of four cells and a small decrease for one cell. In summary, the optimized models show a much closer match to the experimental data as reflected in the individual voltage and current traces as well as in the prediction error. In addition, the optimization identified similar trends in the underlying channel conductance values for different cells from a particular region in the heart. Considered together with the demonstration that the approach accurately identifies model parameters (Figs 2–4), these findings suggest that the approach significantly improves the fidelity of the model for cellular data, relative to the published generic model.
In cardiac modeling, a single action potential or biomarkers derived from it such as amplitude and duration, is often used as a minimal objective for model parameterization. Ionic models can be optimized to fit single action potentials using, e. g. , global search heuristics [22,23], but because the optimization problem is overdetermined, fits may be improved when adding more data, such as data recorded at multiple pacing rates [19,22]. In fact, relative to a single action potential, more complex driving protocols have the potential to dramatically improve parameterization by creating target objectives that are richer in information. On the other hand, to be experimentally feasible, protocols have to be relatively short in duration due to the inevitable current rundown that occurs in patched myocytes, even when using perforated patch. As a compromise, we utilized a stochastic stimulation protocol because it rapidly samples the rate-dependence of the action potential. In addition, irregular excitation patterns are present in many cardiac arrhythmia; thus models tuned to aperiodic excitation patterns are inherently better suited for modeling irregular arrhythmia. In our simulations, we found that the estimates for IKr and IKs were improved the most by the stochastic stimulation objective. During a single action potential, IKr and IKs have similar and compensatory effects (S4 Fig), which impedes estimation of their conductances. In contrast, stochastic stimulation more thoroughly explores their kinetics, thereby revealing small differences throughout the protocol, resulting in a more accurate estimation. Although the stochastic stimulation protocol led to at most a modest improvement in the parameter estimation for the remaining parameters, the prediction error was reduced by an order of magnitude, compared to using a single action potential (Fig 2). Thus, significant model improvement is obtained through the use of a dynamically rich objective, as this helps the optimization avoid the false alternatives that can appear to fit well when dynamically sparse data are used for fitting. In addition to such current-clamp experiments, currents recorded during voltage clamping add additional data to improve fitting and optimization [21,31,32]. While our multi-step voltage clamp protocol alone is very useful for estimating many of the parameter values (Fig 4B), it tends to generate models that fail to predict novel stochastic pacing data well (Fig 4C), which is unsurprising given that it does not train the models according to membrane potential. In our simulations, the addition of the multi-step voltage clamp objective to the stochastic current-clamp stimulation objective enhanced the quality of the parameter estimation compared to using only stochastic stimulation (Fig 4B). This improvement was the result of: (i) some parameters being estimated accurately by the voltage-clamp protocol and (ii) information on two, rather than a single, state variable putting more constraints on the parameter values [36]. In particular, estimates for all nine parameters became centered on their baseline values and the prediction error dropped by another order of magnitude relative to that of stochastic pacing alone. Finally, the iterative optimization approach [31] refined our in silico parameter estimation by decreasing the spread of the returned parameter sets, which caused the prediction error to again decrease by an order of magnitude. Generic models have the advantage that direct comparisons can be made among different simulation studies. However, when comparing a generic model such as the out-of-the-box FR model to our experimental data, there are substantial differences, which likely would cause inaccurate predictions if simulating, e. g. , effects of pharmacological agents or genetic variations. For one, there are clear distinctions in action potential morphology, e. g. , in the plateau phase (Fig 5). This difference in plateau phase most likely explains the method’s downscaling of the IKp conductance. Our recorded action potentials are also of considerably longer duration, which is consistent with the finding of a much reduced IKs in the voltage-clamp experiments. The step to -120 mV in the voltage-clamp protocol induced a much larger current in the experiments than in the FR model and IK1 conductance was increased accordingly in all four cells. These consistent changes in voltage traces and currents between our cells and the FR model may be due to lab-to-lab variability and to the fact that the FR model is not region-specific. Despite such consistent changes, the parameterization also points to important cell-to-cell variability, in particular for the INa conductance, which is increased in three cells and decreased in one. In neuronal modeling, it has become clear that different combinations of conductance parameter sets can give rise to the same activity pattern and that using average values of the conductances may fail to generate that pattern [10,11,37]. The differences in cell-to-cell variation in current densities have been linked to mRNA expression differences or post-translational modifications [38,39]. The extent to which such variation occurs in healthy cardiomyocytes remains to be seen, but some examples of functional coupling among ionic currents in perturbed systems have been described [40–42]. This failure-of-averaging concept may also extend to cardiac tissue: although intrinsic cellular heterogeneity tends to be smoothed out when myocytes are electrically coupled, coupled cells do not necessarily behave like their average. For example, a myocyte with intrinsically shorter action potential duration may promote repolarization in a cell pair [43]. Also, a range of synchronization patterns have been described in coupled pacemaker cells [44]. Thus, there may be important utility to developing cell-specific models. Indeed, cardiac cell-specific models have a range of potential application areas. First, the models can obviously be used to study cell-to-cell variability [45]. Second, in the clinic, inter-subject variability can lead to response differences among patients to pharmaceutical treatment. A dramatic example of this variation is the response to IKr-block, which can vary from minor changes in the electrocardiogram to ventricular tachyarrhythmias [46]. Understanding and predicting this variability is an important step towards patient-specific treatment. In turn, model optimizations such as those developed here represent an advancement towards patient-specific prediction. Finally, multiple models could be grouped into a heterogeneous population and used to generate more realistic responses than those of a randomly-generated population [47,48]. The developed protocols allow accurate estimation of nine conductance/flux parameters. To characterize a single cell more thoroughly, additional flux parameters could be included (e. g. , those describing the sodium/calcium exchanger and the sodium/potassium pump), but as inclusion of more parameters makes the optimization problem harder, this may necessitate tweaking of the methods described here. As detailed below, possible strategies for improvement of the parameter estimations are: 1) improving the stochastic stimulation and voltage-clamp protocols; 2) adding measurements of different state variables during the same protocols (e. g. , intracellular calcium or membrane resistance); 3) incorporating altered solutions and/or ion channel blockers to improve isolation of individual currents [49]; 4) including ion channel kinetic parameters in the optimization; or 5) including relative weights for the current and the voltage contributions to the summed error. Our multi-step voltage-clamp protocol effectively isolates IKs, ICaL, and IK1. An improved voltage-clamp sequence that isolates the remaining currents could improve estimation of their conductance/flux parameters. We designed the voltage clamp steps based on a priori knowledge of the current-voltage (IV) relations in guinea pig ventricular myocytes. As a way to design better protocols, an automated optimization approach may be feasible, i. e. , an optimization of the optimization protocol. Further, differences in structure, channel kinetics and IV-relationships between model and experiment are likely to result in less accurate parameter estimations [31] and may underlie the deviations between fit and experimental data during voltage clamp (Fig 5 and S7–S9 Figs). Adding parameters describing ion channel kinetics to the optimization process would likely improve the fits and predictability, but would almost certainly necessitate longer voltage clamp protocols [6,16]. As channel kinetics are not expected to vary substantially among cells of the same type, a possible strategy is to first parameterize average channel kinetics in a cell population, then apply our method to derive cell specific models. Additionally, improvement could likely be gained by simultaneously recording calcium fluorescence and adding that to the objective function [36], a strategy with merits illustrated by Fig 1 of Ref. [26]. As expected, local sensitivity analysis on simulated calcium traces demonstrates that they are most sensitive to changes in ICaL, IpCa, and JSERCA (S5 Fig), which leads to the speculation that the estimation of these parameters could improve. Inclusion of calcium data may also allow determination of INaCa, which depends on and influences both intracellular calcium and transmembrane potential. Incorporation of membrane resistance in the objective function would also be expected to improve the fitting, as shown in recent work by Kaur et al. [25]. A potential caveat in such multi-objective optimization is that simultaneous good fits are not always achievable, necessitating trade-offs between the different objectives. In that case, balancing which objective (s) to prioritize would be application dependent. Although we allow a generous range for the conductance parameters (0. 01–299% of baseline), some parameters did reach the bounds when fitting the experimental data (Fig 6). Increasing the range will likely require running the GA optimization with a larger population size or for more generations, as will including additional parameters. The main computational cost of the GA is that of simulating the individual models. As this process is inherently parallel, it is straightforward to take advantage of parallel computing. Future implementations could decrease run time by utilizing a GPU, on which optimization for neuronal data has been shown to be feasible [32]. Finally, although the four cells tested in this study provide a strong proof-of-concept for the approach, to further develop the method, it could be applied to a larger number of cells. In the novel approach developed here, cell-specific cardiac models are developed by coupling complex electrophysiology protocols with genetic algorithm parameter fitting. Neither the electrophysiological data (which are too complex to fit by hand), nor the fitting algorithm, would offer much advantage alone. However, merging the two enables markedly improved models that can more accurately simulate dynamically rich cardiac dynamics than can models developed using traditional approaches. Given the widespread use of ionic cardiomyocyte models in investigating arrhythmogenesis, there is utility in models that are better at reproducing such rich electrophysiological dynamics, which are more representative of the complex dynamics that are often inherent to arrhythmias. In addition to improving model fidelity generally, because this approach can be used to generate a model from a single cardiac myocyte, it may be useful in applications ranging from studying the implications of cell-to-cell variability to the prediction of intersubject differences in response to pharmacological treatment.
All animal care and handling for this study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Institutional Animal Care and Use Committee of Weill Cornell Medical College (protocol number: 0701-571A). The cardiac guinea pig model developed by Faber and Rudy (" FR model" ) [34] was used. Model intracellular and extracellular ionic concentrations were set to the values used in our in vitro experiments (see below), after which the model was simulated to a steady state in current clamp mode for 1800 beats at a pacing cycle length of 500 ms. Stimuli were square pulses of 1 ms duration and -40 A/F amplitude. Stochastic stimulation sequences were 5 s in duration. Stimulation times were randomly drawn from a uniform distribution with a range of 100–700 ms. Stimulation times for the optimization sequence were: 216. 48,623. 36,764. 64,1101. 12,1790. 16,2073. 10,2642. 28,3183. 10,3786. 07,3959. 02, and 4579. 72 ms (Fig 2A). Prediction sequence stimulation times were: 247. 40,705. 30,1020. 60,1347. 90,1462. 60,1705. 60,2038. 30,2546. 70,3150. 30,3706. 70,3939. 70,4077. 80, and 4645. 50 ms (S6 Fig). After each parameter change during the GA optimization (details below), current clamp was simulated for 9 beats at a static pacing interval of 500 ms to dampen transients. The voltage clamp protocol was designed on general, a priori, IV relations for the individual channels (e. g. , IK1 is the predominately active current at a holding potential of -120 mV). The 6000 ms protocol was composed of the following steps: 50 ms at -80 mV, 50 ms at -120 mV, 500 ms at -57 mV, 25 ms at -40 mV, 75 ms at +20 mV, 25 ms at -80 mV, 250 ms at +40 mV, 1900 ms at -30 mV, 750 ms at +40 mV, 1725 ms at -30 mV and 650 ms at -80 mV (Fig 3). The contribution of individual membrane currents to the voltage clamp protocol was evaluated using the FR model and the following equation: Contribution=100%⋅∑t=t1t2|Ix (t) |∑j=1N∑t=t1t2|Ij (t) | (1) Eq 1 calculates the percentage contribution of the absolute individual current (Ix) relative to the absolute sum of all N currents for all time points during one of the holding potentials (t1 to t2). This calculation was done for all model currents at all holding potentials. During GA optimization, the simulated multi-step voltage clamp is preceded by 5 s holding at -80 mV to allow the model to settle after parameter changes. Multiple global search heuristics have been applied to electrophysiology models, including gradient-based descent [13,19–21], simulated annealing [50], particle swarms [18,24], and genetic algorithms [22,23,25]. We chose a genetic algorithm as it is effective for a range of the number of parameters [50], is computationally simple and readily parallelizable, and has been shown to be successful at optimizing sophisticated ionic models to experimental data [22,23,25]. The GA used in this study was originally developed by Sastry [33] and was used with the settings for selection, crossover, mutation, and elitism strategy as in Bot et al. [23]. Compared to the study of Bot et al. , we increased the parameter search range to 0. 01–299% of the baseline model values. This larger range, in combination with the increased number of parameters and the diminished requirements for computation speed relative to the Bot et al. study, caused us to enlarge the population size to 500 and raise the number of generations to 100, based on test runs showing consistent convergence when using these values. Because the GA is inherently stochastic, it was run 10 times per optimization problem. In addition, an iterative approach was implemented, based on the study of Hobbs and Hooper [31]. For each parameter, the iterative approach uses the span of the 10 best individuals from the first 10 GA runs as the search boundaries for a second set of GA runs. With these new search ranges, the GA was again run 10 times with a population of 500 individuals for 100 generations. We used mean squared differences for the objective functions (errors) that the GA works to minimize: E1=∑t=tIC, starttIC, end (Vtarget (t) −Vindividual (t) ) 2 (2) E2=∑t=tIC, starttIC, end (Vtarget (t) −Vindividual (t) ) 2+∑t=tVC, starttVC, end (Itarget (t) −Iindividual (t) ) 2 (3) where E1 is the objective function when only current clamp data (i. e. , stochastic stimulation or a single action potential) was fit, and E2 is the objective function for the combined stochastic stimulation and multi-step voltage clamp protocol. In both equations, Vtarget is the membrane potential during current clamp of the target (i. e. , either the simulated nominal model or the experimental data), and Vindividual is the membrane potential of a simulated individual. Itarget and Iindividual are the current responses during voltage clamp of the target and a simulated individual, respectively. Errors are summed over the entire duration of the protocols. Although of different units, the voltage clamp and the current clamp components to E2 were simply summed into a single objective (Eq 3) as we expect them to be minimal for the same range of parameters, rather than being competitive as in typical multi-objective optimization. E2 is therefore unitless. The estimated model parameters are the maximal conductances of INa, ICaL, ICaT, IK1, IKr, IKs, IKp, and IpCa, and the maximal flux of JSERCA. Optimizations were run on a 3. 2Ghz Intel Xeon W3670 6-core, 6GB memory, machine and took approximately 8 hours per run for the iterative approach and then combined stochastic pacing and voltage clamp protocol. Two sample t-tests were performed with a significance level of 0. 05. Numbers and error bars indicate average ± standard deviation. Guinea pigs (n = 4) were anesthetized using an intraperitoneal injection with Euthasol (Virbac Corporation, Fort Worth, TX), 550 mg/kg. Excised hearts were then Langendorff retrograde perfused, and myocytes were isolated from the base (top 1/3) of the left ventricle through enzymatic digestion. Myocytes were stored in Dulbecco’s Modified Eagle Medium (DMEM) with 5% fetal bovine serum (FBS). Amphotericin-B (Sigma-Aldrich Corp. , St. Louis, MO; 480 μg per 1 ml pipette solution) perforated patch clamp technique was used to record cellular action potentials. Bath solution contained (in mmol/l) 139. 4 NaCl, 5. 4 KCl, 1. 0 MgSO4,10. 0 Hepes, 10. 0 dextrose, 2. 0 CaCl2, pH 7. 35 with NaOH, osmolality 310 ± 3 mmol/kg. Intracellular solution contained (in mmol/l unless otherwise noted) 125 KCl, 10 NaCl, 5. 5 dextrose, 0. 5 MgCl2,11 KOH, 10 Hepes, 10 μmol/l CaCl2, pH 7. 1 with HCl, osmolality 295 ± 3 mmol/kg. Recordings were performed at 35°C. Patch-clamp measurements were recorded using an Axopatch 200A amplifier (Molecular Devices, Sunnyvale, CA). The Real-Time eXperiment Interface [RTXI; rtxi. org; [51,52]] software platform was used to control the amplifier and record data. Cells were initially paced in current clamp mode at a BCL of 500 ms to steady state (500–1000) beats using suprathreshold square pulses of 1 ms duration. Next, the optimization and prediction stochastic stimulation sequences were applied. Amplifier mode was then switched to voltage clamp and series resistance measured (4–8 MΩ) and compensated for (70–90%). The multi-step voltage clamp protocol was then applied in triplet. Holding potentials were corrected for a liquid junction potential of -3 mV. The magnitude of the Donnan equilibrium was estimated to 0 mV using the IV-curve of INa and therefore not corrected for. To remove stimulus artifacts from current-clamp traces, data from a 1. 3 ms window following the start of each stimulus were excluded from the optimization. In addition, voltage-clamp data from a 1. 2 ms window following each potential change were excluded from the GA optimization because of the capacitance transient. From the set of three voltage-clamp trials, the current response trace with the shortest time to peak INa (step to -40 mV at 600 ms) was selected for each cell.
Mathematical models of cardiac cell electrophysiology are widely used as predictive and illuminatory tools, but have been developed for decades using a suboptimal process. The models are typically constructed by manual adjustment of parameters to fit simple data and therefore often underperform when used to predict complex behavior such as arrhythmias. We present a novel method of model parameterization using automated optimization and dynamically rich fitting data and then demonstrate that this approach is better at finding the “real” model of a cell. Application of the method to cardiac myocytes leads to cell-specific models, which may enable well-controlled studies of both cellular- and subject-level population heterogeneity in disease propensity and response to therapies.
Cell-Specific Cardiac Electrophysiology Models