Cell-Specific Cardiac Electrophysiology Models

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.


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][3][4][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: 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 data [9].
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.
3. 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-tocell 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].
4. 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.
5. 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][13][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.
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][16][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][19][20] and more physiologically detailed ionic models [21][22][23][24][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][28][29][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.

GA optimization using a single action potential
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 (I Na ), the L-type calcium current (I CaL ), the T-type calcium current (I CaT ), the inwardly rectifying potassium (I K1 ), the rapid and slow delayed rectifier potassium currents (I Kr and I Ks ), the plateau potassium current (I Kp ), and the sarcolemmal calcium pump current (I pCa )] and the maximal flux of the sarcoplasmic reticulum Ca 2+ -ATPase (J SERCA ) 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.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 ).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
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 I Kr , I CaL , and I Ks (Fig 2C).However, as shown in Fig 2C, a few parameters remain incorrectly estimated (maximal conductances of I CaL and I Ks ), and some are estimated with a large spread (maximal conductance values of I Kp , I pCa , I CaT , and I Kr ).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) 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.
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 I K1 , I CaL , and I Ks 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 Indeed, using the voltage clamp protocol as the objective during an optimization recovers the conductances for I K1 and I Ks very accurately (Fig 4B).The estimation of the I CaL conductance is very close to 1, but is slightly overestimated in all runs.Less predictively, I pCa and I Kp were also estimated more precisely than during stochastic pacing alone (Fig 4B).However, a few conductances were estimated poorly (in particular J SERCA and I CaT ) 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 (I Ks , I K1 , and I pCa ) 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 I Na and voltage clamp for I Ks ).However, for other currents, in which both individual protocol segments result in off-target outcomes, the combined protocol produces estimates spanning 1 (e.g., I Kp ).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 I Na .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  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.
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).

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
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 ratedependence 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 I Kr and I Ks were improved the most by the stochastic stimulation objective.During a single 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 currentclamp 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
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 I Kp conductance.Our recorded action potentials are also of considerably longer duration, which is consistent with the finding of a much reduced I Ks 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 I K1 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
Despite such consistent changes, the parameterization also points to important cell-to-cell variability, in particular for the I Na 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][41][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, intersubject variability can lead to response differences among patients to pharmaceutical treatment.A dramatic example of this variation is the response to I Kr -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 patientspecific 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
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 I Ks , I CaL , and I K1 .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 I CaL , I pCa , and J SERCA (S5 Fig), which leads to the speculation that the estimation of these parameters could improve.Inclusion of calcium data may also allow determination of I NaCa , 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.

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

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.

Multi-step voltage clamp protocol
The voltage clamp protocol was designed on general, a priori, IV relations for the individual channels (e.g., I K1 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: Eq 1 calculates the percentage contribution of the absolute individual current (I x ) relative to the absolute sum of all N currents for all time points during one of the holding potentials (t 1 to t 2 ).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 gradient-based descent [13,[19][20][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: Although of different units, the voltage clamp and the current clamp components to E 2 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.E 2 is therefore unitless.
The estimated model parameters are the maximal conductances of I Na , I CaL , I CaT , I K1 , I Kr , I Ks , I Kp , and I pCa , and the maximal flux of J SERCA .
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 ± standard deviation.

Electrophysiology experiments
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.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 MO) 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 I Na and therefore not corrected for.

Post-processing 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 I Na (step to -40 mV at 600 ms) was selected for each cell.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 are directly calcium-related (I CaL , J SERCA , and I pCa ).

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

Fig 1 .
Fig 1. Progression of GA estimation.(A-C) The GA is initialized with 500 random individuals, i.e., model instantiations, in generation 0. Each individual is simulated according to the protocol, here a single action potential, and its error is calculated (sum of squared errors between model output and FR model target objective; Eq 2 in the Methods).Left columns show action potentials generated by three different generation 0 model instantations (traces are colored according to their error and color bar in panel F) compared to the baseline FR model (black).Right column bar graphs indicate the scaling of the nine model parameters for each individual, with a scaling of 1 representing the original FR model value.Parameters 1 through 9 correspond to: I Na , I CaL , I CaT , I K1 , I Kr , I Ks , I Kp , I pCa , and J SERCA , respectively.(D-G) With progression through the generations, individual action potentials become more similar to the optimization objective and errors decrease accordingly.At generation 100, the overall best individual and the FR model appear very similar, although the bar graph indicates differences among the parameters (E).doi:10.1371/journal.pcbi.1004242.g001 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.

Fig 2 .
Fig 2. Stochastic stimulation improves parameter estimation and predictability.(A) The stochastic stimulation optimization sequence (black) and the best fit from 10 GA runs applied to fit it (light blue dashed line) show the quality of the GA fit.(B) Zoom-in of shaded region in (A) shows fitness of closely coupled action potentials.(C) Using the stochastic stimulation optimization sequence in the objective function causes some improvement in parameter recovery (light blue circles) compared to using a single action potential (green squares), most notably for I Kr , I CaL , and I Ks conductances.Symbols indicate the best individual from each of the 10 individual GA runs, error bars give the mean (black) ± standard deviation (gray), and the dashed line at parameter scaling 1 indicates the original FR model parameter values.(D) When subjected to a novel stimulation sequence, the error between the model fit based on the stochastic stimulation sequence and the FR model (light blue circles) is significantly smaller (p-value of 0.0001) than the error between FR model and the best fit based on a single action potential only (green squares), i.e., the stochastic stimulation objective leads to models with higher prediction accuracy.doi:10.1371/journal.pcbi.1004242.g002

Fig 3 .
Fig 3. Multi-step voltage clamp protocol.(A) Voltage clamp protocol (blue) and current response (red).(B) Close-up of shaded region in (A) where voltage steps change rapidly.(C-G) Percentage contributions of the individual currents during particular parts of the protocol isolates I K1 (C), I CaL (E) and I Ks (F and G) well.doi:10.1371/journal.pcbi.1004242.g003

Fig 4 .
Fig 4. Improvement of the optimization results with inclusion of voltage-clamp protocol and use of iterative optimization technique.(A) Error reduction during optimization with combined stochastic pacing and voltage clamp objective.Circles indicate contribution to the summed error from each of the 500 individuals during stochastic pacing (blue) and voltage clamp (red), while the correspondingly colored traces give the population means.The summed error of the best individual is shown as the black trace.(B) Improvements in the parameter estimation from using stochastic stimulation (light blue circles) or voltage clamping (green squares), to using the combined stochastic pacing and voltage clamp protocol (orange diamonds), and to using the iterative optimization (magenta triangles) as visualized with estimation results centered around 1 and tighter error bars.Symbols indicate results from the individual runs; error bars give the mean (black) ± standard deviation (gray).(C)The prediction error is large for the voltage clamp protocol alone (green squares), which does not train models according to membrane potential.Adding the voltage clamp protocol to the stochastic pacing protocol (i.e., combined protocol; orange diamonds) gives better predictions compared to stochastic stimulation alone (light blue circles; p = 0.00005).The prediction error is further reduced when using the iterative optimization approach (magenta triangles; p = 0.0003).doi:10.1371/journal.pcbi.1004242.g004

Fig 5 .
Fig 5. GA voltage and current recovery for in vitro myocyte compared to FR model.(A) Stochastic stimulation optimization sequence shows the difference between the experimental voltage response (red trace) and the original FR model (black trace).The GA fitted model (blue, dashed trace) matches the experimental action potentials closely.(B) The multi-step voltage clamp protocol reveals that the experimental current data is also very different from the original FR model, in particular during the first I K1 -inducing step which results in a much larger current in the experimental recording, and during I Ks -inducing steps which trigger much larger currents in the model.Again, the experimental data is well fit by the optimized model.(C) The optimized model also predicts the response to a novel stochastic stimulation sequence well.This is reflected in the prediction error (D), which is reduced for the optimized models for all four cells (colored circles; standard deviation of the 10 best individuals for each cell falls with the circles) compared to the FR model (black circles).The data shown in (A-C) is from Cell 2. In (A) and (B) stimulus artifacts and capacitative currents were removed (see Methods), but data sets were plotted as continuous traces to ease visualization.doi:10.1371/journal.pcbi.1004242.g005
(PDF) S6 Fig. Stochastic pacing prediction.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.(PDF) S7 Fig. Experimental data fit, cell 1.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 Fig 5), but data sets were plotted as continuous traces to ease visualization.(PNG) S8 Fig. Experimental data fit, cell 3. Stochastic pacing and voltage clamp fits (blue) of the experimental data of cell 3 (red) compared to the original FR model (black).(PNG) S9 Fig. Experimental data fit, cell 4. Stochastic pacing and voltage clamp fits (blue) of the experimental data of cell 4 (red) compared to the original FR model (black).(PNG) is the objective function when only current clamp data (i.e., stochastic stimulation or a single action potential) was fit, and E 2 is the objective function for the combined stochastic stimulation and multi-step voltage clamp protocol.In both equations, V target is the membrane potential during current clamp of the target (i.e., either the simulated nominal model or the experimental data), and V individual is the membrane potential of a simulated individual.I target and I individual 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.