Abstract

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.

Author Summary

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.

Introduction

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].

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: 1. Because individual quantification of all membrane currents requires many experiments, model-development data are typically taken from multiple laboratories, often using different experimental protocols with varying conditions such as temperature and solution composition that directly influence the ionic currents [7,8]. Such variations may be taken into account by the modeler, but extrapolations to other conditions are often based on sporadic 2. Data may originate from cells from different locations within the heart, or even different species, Which can be problematic because of regional/ species heterogeneities in electrophysiological properties [8]. The data may also come from animals or patients With divergent characteristics, such as gender and age, that influence electrical activity.

An additional cause of inconsistency comes from averaging over multiple experiments of both the voltage clamp data used to build the model and the functional data used to parameterize it. Even in myocytes from the same region of the heart, there is considerable cell-to-cell variability in the action potential, presumably stemming from different levels of ion channel densities. Thus, models create the ‘average’ behavior for a cardiac myocyte, although no single cell necessarily behaves like the average model [10,11].

Another problematic issue arises when the data and equations for the individual currents are combined to build a composite model. This step typically requires tuning of model parameters to reproduce Whole-cell behavior; this tuning is usually done manually in a laborious, iterative tweaking process, Which ends When the model output (e.g., an action potential) subjectively is deemed to adequately match the experimental counterparts. Such manual adjustments will almost certainly not result in the best possible fit to the data, as the approach can explore only a small fraction of parameter space.

Given the relatively simple dynamics (e.g., often a single action potential) to which the models are tuned, models often fare poorly When attempting to simulate more complex behaviors (e.g., fibrillatory dynamics, drug block, or pacing interventions) [7,12—14]. This is unfortunate as one of the key goals of cardiac modeling is to predict and give mechanistic insights into arrhythmias, Which often occur at excitation patterns different from those used to construct the models.

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.

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.

Results

GA optimization using a single action potential

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 (1K1), the rapid and slow delayed rectifier potassium currents (1Kr and IKS), the plateau potassium current (IKp), and the sarcolemmal calcium pump current (IpCa)] and the maximal flux of the sarcoplasmic reticulum Ca2+-ATPase (ISERCA) were estimated using a GA technique.

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].

Stochastic stimulation protocol improves model fit and predictability over single action potential protocol

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).

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 1K3 (Fig 2C). However, as shown in Fig 2C, a few parameters remain incorrectly estimated (maximal conductances of ICaL and 1K5), 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 deter-mined—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) (SI—$4 Figs and 81 Text).

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.

Extending the protocol: Adding multi-step voltage clamp data

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.

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 1K1, ICaL, and 1KS 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.

The combined stochastic current and multi-step voltage clamp protocol improves parameter estimation

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 ISERCA 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).

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.

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, 1K1, and IpCa) did not show improvement with the combined protocol.

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.

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.

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.

Dynamic electrophysiology protocols and optimization improve model fit to in vitro experimental 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.

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).

Parameter estimation shows changes compared to FR model and variability among individual cells

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 ISERCA is slightly reduced. In contrast, for all four cells, maximal conductance of IKr and 1K1 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 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.

Discussion

To overcome the limitations inherent to traditional cardiac model construction (most notably manual parameter adjustment and use of averaged, dynamically limited data), we developed a novel approach for parameter estimation that combines an electrophysiology protocol that is rich in dynamic information, short in duration, experimentally feasible, and does not require the use of drugs or special solutions, with a parallel computational parameter tuning algorithm (GA). The protocol was first tested computationally, which showed reliable parameter estimation. We then applied the protocol in vitro to guinea pig basal left ventricular myocytes. Compared to the baseline guinea pig FR model, the optimized cell-specific models showed a significantly improved fit to the experimental data. Furthermore, because our model enables validation on data from the same cell for which a model was optimized, we were able to demonstrate that the cell-specific models are markedly better at predicting the response to novel stimulation sequences than was the generic model.

Complex driving protocols and objectives

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, f1ts 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.

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.

Improvement in model parameterization for intact cardiac myocytes

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 methods 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 1KS 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 1K1 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.

Cell-specific models

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.

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] .

Limitations and potential improvements

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.

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.

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.

[26]. As expected, local sensitivity analysis on simulated calcium traces demonstrates that they are most sensitive to changes in ICaL, IpCa, and ISERCA (SS 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 tradeoffs between the different objectives. In that case, balancing which obj ective(s) to prioritize would be application dependent.

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.

Conclusions

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 myo-cyte, 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.

Methods

Ethics statement 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).

Computational modeling

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 protocol

500 ms to dampen transients.

Multi-step voltage clamp protocol

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: 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.

Genetic algorithm optimization

Multiple global search heuristics have been applied to electrophysiology models, including gra-dient-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].

[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: 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. 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.

Statistics

Two sample t-tests were performed with a significance level of 0.05. Numbers and error bars indicate average i standard deviation.

Electrophysiology experiments

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).

Cells were initially paced in current clamp mode at a BCL of 500 ms to steady state (500—1000) beats using suprathres-hold 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 M9) 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.

Postprocessing of experimental data

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.

Supporting Information

Supplemental results. Methods and results of sensitivity and correlation analyses. (PDF)

Parameters were scaled to 80, 90, 95, 105, 110 and 120% of their published values and the sum of squared errors was calculated and Visualized here as the sensitivity. For each parameter, the effect of the scaling is given from small to large parameter scaling, i.e., from 80—120%. See 81 Text for details. (PDF)

Multi-parameter sensitivity analysis was carried out as detailed in 81 Text. Blue symbols indicate parameters with a statistically significant effect on the model output While the red symbols indicate the parameters with a nonsignificant effect on the model output (ICaT, 11(1), and IpCa). The dashed line indicates the largest sensitivity of the nonsignificant parameters and is Visualized as a threshold value for sensitivity. (PDF)

The figures give the squared coefficient of variation (standard deviation normalized to the mean) for each conductance/flux parameter during the optimization process, averaged for the 10 GA runs.

See 81 Text for details. (PDF)

Colors represent the value of the correlation between two parameters. Symbols indicate statistical significance. See 81 Text for details. (PDF)

Parameters were scaled to 80, 90, 95, 105, 110 and 120% of the published value and the sum of squared errors (using intracellular calcium concentration rather than transmembrane potential or total current in Eqs 2 and 3) was calculated and Visualized here as the sensitivity. For each parameter, the effect of the scaling is given from small to large parameter scaling, i.e., from 80—120%. The calcium signal is most sensitive to parameters that

A) Prediction sequence used to calculate prediction error. Best individual from stochastic pacing GA optimization runs (blue, dashed) and FR model (black) show close correspondence. B) Prediction error calculation for the best individual from 10 GA optimization runs using a single action potential (green) or stochastic pacing (blue). FR model simulation (objective) is given in black. The individual from the stochastic pacing runs matches the FR objective more closely.

Stochastic pacing and voltage clamp fits of the experimental data of cell 1. The figure shows the best individual from 10 GA runs using the iterative approach (blue), the original FR model (black) and the experimental data (red). The GA fit shows a closer match with the experimental data than the FR model. Stimulus artifacts and capacitative currents were removed (as in Egg), but data sets were plotted as continuous traces to ease visualization.

Stochastic pacing and voltage clamp fits (blue) of the experimental data of cell 3 (red) compared to the original FR model (black). (PNG)

Stochastic pacing and voltage clamp fits (blue) of the experimental data of cell 4 (red) compared to the original FR model (black).

Acknowledgments

The authors thank Erica Bishop and Lala Tanmoy Das for help with the isolation procedure and preparation of the solutions.

Author Contributions