Efficient parameterization of cardiac action potential models using a genetic algorithm

Many models of cardiac action potentials exist, but the large level of biological variability stemming from differences including species, region of the heart, age, sex, general level of health, and natural intra- and inter-subject variability, combined with differences in conditions between laboratories, makes it extremely difficult to match a model as published to specific experimental data. Nevertheless, there is a growing need for quantitatively accurate models that can be used reliably in applications from personalized health care^8–16 to drug development and regulatory approval.¹⁷ Thus, it is essential to develop tools to find parameter values that allow a model to match specific experimental data.

Genetic algorithms (GAs) have the potential to help with parameterization of these models because they are designed to explore large portions of parameter space and to reduce the likelihood of being caught in a local optimum. Inspired by strategies used in evolution, GAs utilize a population of candidate solutions that evolve over multiple generations through the processes of selection, crossover, and mutation. The main idea is to encourage improvement of good solutions while maintaining sufficient diversity to consider other potentially good solutions that may be located farther away in parameter space. GAs have been used in the past to find parameter sets for cardiac models that reproduce experimental data. Syed et al.¹⁸ used a genetic algorithm to find values for nine ion channel conductance parameters of the human atrial model of Nygren et al.¹⁹ to fit action potentials obtained from the human atrial model of Courtemanche et al.,²⁰ which has a different action potential morphology at long cycle lengths (CLs), as well as experimental human atrial potential recordings, where data from multiple pacing frequencies were used. Kherlopian et al.²¹ also used a genetic algorithm to fit five conductances of the Luo-Rudy I model²² and eight conductances of the model of Nygren et al. to guinea pig ventricular action potentials. Bot et al.²³ applied the work of Kherlopian et al. to fit six conductance parameters of a neonatal mouse ventricular^24,25 model to match experimental data. Later, Groenendaal et al.²⁶ extended this work to fit guinea pig ventricular myocytes with additional data from voltage-clamp experiments. In all of these cases, however, only as many as nine parameters were fit, representing a small fraction of the full complement of model parameters.

Here, we systematically study an alternative approach in which we use a genetic algorithm to fit all parameter values of the models of Fenton-Karma²⁷ and Bueno-Orovio et al.,⁵ which have been shown previously to be flexible in reproducing the dynamics of other models^5,27,28 as well as experimental and clinical data^{5,8,15,27,29–33} using only three or four state variables. Such models are not fully mechanistic in that they do not reproduce the detailed behavior of ion channels, pumps, and exchangers, but they are quite valuable for large-scale tissue simulations. We allow our algorithm to determine values for all 13 parameters of the three-variable (3 V) Fenton-Karma model and all 27 variables for the four-variable (4 V) model of Bueno-Orovio et al. (we only impose that the value corresponding to resting potential to be 0, which is the value used unless it is important to model spatially varying differences in resting membrane potential³⁴). We use a fitness function that prioritizes both reducing curve error and matching action potential duration.

We find that a genetic algorithm performs well for fitting both model-derived and experimental data in many cases, although dynamics can vary significantly for CLs not included during the fitting process, especially during alternans cases. In addition, we find that there can be a good deal of variability in the specific parameter values chosen, which points to underlying issues of parameter identifiability and sensitivity for these models.

A genetic algorithm (GA) is a biologically inspired method for solving optimization problems. It considers a pool of candidate solutions (parameter sets), called a population, and allows it to evolve over some number of generations, or iterations. Each new generation is created by repeatedly applying the basic operators GAs use to explore their search space: selection, crossover, and mutation. In each generation, every member of the population is assigned a fitness value using an objective function that quantifies how well a given parameterization solves the problem. The fitness function is one of the most important parts of the algorithm, and its specifications influence the overall performance at every level. More information regarding these algorithm components is available in the supplementary material.

The following default algorithm settings are used unless noted otherwise: tournament selection with tournament size of 5; one-point, one-child crossover; 2500 population members; 10 generations; Gaussian mutation with a rate of 0.1%, mean mutation size of 5% of the bounds for that parameter, and standard deviation as 1/3 of the mean mutation size; and elitism (see supplementary material). Bounds for parameters are set to be quite broad by including the ranges of previously published parameterizations by Fenton and Karma²⁷ and Fenton et al.³⁵ for the 3 V model and by Bueno-Orovio et al.⁵ for the 4 V model. The initial population is generated by selecting values within the established bounds for each parameter according to a uniform distribution. Using an approach to guarantee parameter space coverage like Latin hypercube sampling did not show any significant difference in results and slightly extended runtime.

To understand the performance of the genetic algorithm in different contexts, we use two flexible models of cardiac cellular electrophysiology: the three-variable (3 V) model of Fenton-Karma and the four-variable (4 V) model of Bueno-Orovio et al. These models have been shown previously^{5,27–29,32,35} to be highly useful for reproducing, through different parameterizations, a broad range of dynamics of cardiac cells and tissue obtained experimentally or from other more complex models. Both models include fast inward, slow inward, and slow outward currents roughly corresponding to summary sodium, calcium, and potassium currents.

The 3 V model includes a voltage variable $u$ along with inactivation gating variables $v$ and $w$ for the fast inward and slow inward currents, respectively. Its 13 model parameters can be tuned to reproduce changes in action potential duration as a function of CL, and many of its parameters have specific interpretations within the context of the model (e.g., $uc$ is the threshold for excitation, $ucsi$ is the threshold for activation of the slow inward current, $τd$ represents excitability and governs the upstroke steepness, $τr$ is the time constant governing early repolarization, $τo$ is the time constant governing late repolarization, $τv1−$ is proportional to the minimum diastolic interval, etc.). The model is “semi-normalized” in that the upstroke, but not the full action potential, is normalized to lie between 0 and 1 for long CLs; the voltage $u$ may increase further later in the action potential. A limitation of the 3 V model is that the range of action potential shapes that can be achieved is fairly limited.

The 4 V model extends the 3 V model to improve the realism of action potential morphology by introducing a fourth variable, $s$⁠, which serves as an activation gating variable for the slow inward current. Thus, there is a delay in the development of the slow inward current following the upstroke, which is more faithful to the underlying physiology and can result in a greater variety of more realistic action potential shapes. Action potentials typically are scaled to fall between 0 and around 1.4–1.5. Along with the addition of a fourth variable, the 4 V model allows greater flexibility by roughly doubling the number of model parameters. Tuning its 27 model parameters can be challenging; a gradient-based optimization routine was used to obtain the fittings given in the original publication.⁵ 

We use the explicit Euler method for time integration with a time step of 0.05 ms. Initial values are set to $u=0$ and $v=w=1$⁠; for the 4 V model, we also set $s=0$⁠. To improve code performance, we use a lookup table to calculate values of the model's tanh functions. Specifically, we construct a table of tanh values for arguments ranging from 0 to 4 in increments of 10⁻⁴. Arguments greater than 4 are assumed to have function values of 1. The functions of negative arguments are found as the negations of the functions of their absolute values (because tanh is an odd function). Performance is further enhanced by computing fitness values for the population in parallel using the POSIX threading library.

Both model- and experiment-derived data are used for obtaining parameterizations. For model-generated data, only the output for the voltage variable every 1 ms is used for fitting. For each CL used, all model variables are set to the initial conditions and paced for 10 beats, which is usually sufficient for these models because of the minimal amount of memory present. The last two action potentials (beats 9 and 10) are used for fitting to allow alternans properties to be considered as a part of the fitting process.

Experimental data were obtained from microelectrode recordings of canine cardiac tissue at physiological temperatures. Recordings of the voltage were made at 1 kHz over a downsweep pacing protocol; the tissue was paced to a steady state at each CL included in the protocol. Representative successive action potentials obtained during steady state were extracted for each CL. The data were processed using custom software to remove any drift and stimulus artifacts and to normalize the voltage between 0 and 1. No time-averaging was used in order to preserve the action potential upstroke features; thus, a small amount of noise was present in the signal. To construct bifurcation plots from experimental data, APDs were measured at 90% of repolarization and were averaged over the last 10 beats when alternans was not present and over 10 long and 10 short beats when alternans was present.

APDs are measured in all cases using a threshold value of 90% repolarization, corresponding to $u=0.1$ and $u=0.145$ for the 3 and 4 V models, respectively. Experimental data are rescaled to a maximum of 1.45 for the 4 V model for the longest CL because of its higher plateau value prior to APD measurement. APDs measured within the genetic algorithm implementation are interpolated linearly. CLs resulting in 2:1 block are excluded from the bifurcation plots.

To keep computational costs low, our genetic algorithm uses a single-cell formulation. However, it is known that many phenomena in single cells do not translate straightforwardly to tissue.³⁶ An important difference is in how cells are stimulated: in single cells, a stimulus must be provided directly to the cell using a time-dependent current, often a square pulse lasting 1–2 ms with a magnitude set to twice the minimum value needed to produce a stimulus. In tissue, although some cells may be stimulated directly in the same manner as for isolated cells, most cells are stimulated via electrotonic currents arising from wave propagation through diffusive coupling. The electrotonic current seen by a cell in tissue as it becomes excited differs in several ways from a short-duration square pulse: (i) the current is nonzero for a longer period of time; (ii) the current is not square but rather has a magnitude that changes continuously; and (iii) the current is biphasic, with coupling current changing sign after the cell has become excited and its neighboring cells draw current to become excited themselves.

To improve the likelihood that single-cell simulation results correlate well with findings obtained in tissue, thereby allowing us to keep the algorithm runtimes short, we use a biphasic stimulus current for each model saved from a one-dimensional simulation. Additional information about the biphasic stimulus protocol, including how it can screen out some 0d cases where propagation would fail in tissue, can be found in the supplementary material.

To determine whether two series of algorithm runs used to fit two different data sets resulted in statistically different values for a given parameter, we used a Wilcoxon Rank Sum test with a level of significance of $p=0.01$⁠. The test was applied for all parameters to identify which parameters were statistically significantly different for fittings to the two experimental data sets.

We first study how the algorithm works by considering several cases where only a small number of parameters are fit, which allows us to visualize what happens as the algorithm progresses. We then assess the algorithm's performance for so-called “model recovery” cases where the data used for fitting are derived directly from the model that is fit. Finally, we use the algorithm to obtain parameter fittings for experimental data obtained from microelectrode recordings in cardiac tissue.

When fitting only two or three parameters, it is possible to obtain a useful visualization of what is happening as the algorithm runs. In particular, we can calculate the fitness function used in the algorithm over the entire parameter space, with the bounds of each parameter used in fitting serving as the bounds in the fitness function calculation. Figure 1 shows results for fitting two parameters for each model: $τd$ and $uc$ for the 3 V model and $τsi$ and $τw+$ for the 4 V model. Note that the results shown in these figures are from a single run of the algorithm; because the algorithm includes randomness in producing the initial population as well as in selection, crossover, and mutation, additional runs can give different results. The lowest values of the fitness function generally occur for cases where the parameterization does not result in the production of action potentials or when action potential block occurs. For the 3 V model case, there is a broad range of parameter value combinations that lead to reasonable action potentials and thus relatively high values of the fitness function (over 0.5), which can be seen as the large amount of gray present in the grayscale fitness background. For the 4 V model case, the function is more restrictive, and only a narrow region within parameter space has high values of fitness. In both cases, population members consist of pairs of parameter values and are superimposed as colored dots, with the correct parameter values shown with a black cross. As data to fit, the last two action potentials from a train of ten at CLs of 500, 400, and 350 ms are used using sets 3 V-1 and 4 V-1 (full parameter values are given in the supplementary material).

The initial population at the beginning of generation 1 is fairly evenly distributed. The calculated fitness values are used to choose parents for producing the next generation, and, as expected, the majority of population members selected (magenta) have much higher fitness values than those not selected (blue). Nevertheless, some population members without high fitness values are selected. This behavior is important because while the algorithm should be improving on good fitness values, it also should retain enough diversity to consider alternatives that can allow exploration of other areas of parameter space. In these examples, however, the fitness functions are relatively simple; thus, the algorithm proceeds fairly smoothly toward identifying the global maximum.

Each generation behaves similarly, but with an increasing focus on the area of the fitness function corresponding to the correct parameter values, which have the maximum fitness. Note that the minimum fitness value of two parents is not the minimum of the fitness value for a child of those two parents: thus, it is possible for crossover to produce more diverse individuals than their parents. This feature ensures that more of parameter space can be investigated rather than ruling out regions too quickly.

The action potentials produced by the best parameterizations obtained for each case (not shown) are so close to the target action potentials that it is nearly impossible to distinguish the two. For the 3 V case, the difference between the target value of $u$ and the value of $u$ produced by the parameterization is less than $2×10−5$ for all times; for the 4 V case, the maximum difference is $1.04×10−3$⁠. Differences outside the action potential upstroke and repolarization are negligible.

We also show a case in which three parameters are fit within the 4 V model. Now instead of viewing the fitness function in the three-dimensional parameter space, we simply display several contours using two different views, as shown in Fig. 2. Again, the population initially is spread fairly uniformly, but because of the simplicity of the fitness function, the algorithm quickly converges to the global maximum within the innermost contour. As with the two-parameter cases, the action potentials obtained using the parameterization found are very close to the target action potentials and thus are not shown. Here, the maximum difference in $u$ at any time is $3.11×10−3$⁠.

We note that the fitness functions for smaller number of parameters, such as those depicted in Figs. 1 and 2, tend to be relatively simple. When more parameters are included, the fitness functions can easily become more complicated and may have multiple local maxima that may not be close to each other. Although the visualization approach we have used here cannot be generalized to a higher-dimensional parameter space, similar principles apply.

To test the algorithm more fully, we again consider cases in which we use data sets derived from the model whose parameter values were sought, but this time with all parameters included in the fitting. Thus, 13 parameters for the 3 V model and 27 parameters for the 4 V model are fit (as mentioned earlier, the 4 V model resting membrane potential parameter $uo$ is not fit but is set to 0).

Figure 3 shows the results of fitting the 3 V model to the Beeler-Reuter parameter set (set 3 V-1 in the Supplementary material) for a single run of the algorithm. Action potentials from CLs of 500, 400, and 350 ms are fit. As can be seen, the fitting obtained produces very similar action potentials, with the maximum difference in APD less than 2 ms for all fit CLs. The differences are more noticeable here than when only two model parameters are fit because of the greater complexity of parameter space (which is now 13-dimensional). Nevertheless, good agreement is obtained, even for many unfit CLs. For CLs below the bifurcation point, however, the model is very sensitive to the specific parameterization. Here, disagreement in the bifurcation plot is evident for CLs below 350 ms. In this case, if capturing the alternans dynamics were important, it would be necessary to include data from more small CLs in the fitting process.

The results just described are for a single run of the algorithm. However, the random components of the algorithm cause variations in performance. Viewing results from multiple runs together can give an idea of the consistency of the algorithm in selecting the values of each particular parameter. Figure 4 shows results from 50 runs of the algorithm, specifically as box-and-whisker plots of the parameter values obtained. The algorithm exhibits strong preferences for the values of some parameters. For example, $τr$ is allowed to vary over a very broad range, but it never chooses a value outside the bottom quarter of its range. Similarly, the values for $τw−$ are very well clustered, and $ucsi$ never takes on a value below the top half of its range. The most variability can be seen for the values of parameters $τv1−$ and $τv2−$⁠. This can be explained by the fact that these parameters primarily control the properties of conduction velocity restitution curve slope and minimum diastolic interval, neither of which can be recovered using the data provided. The value of $τw+$ also is quite variable and, as it is related to the maximum action potential duration possible, also would be impossible to deduce from the data used for fitting.

To provide an idea of the dynamical variability of the parameter sets recovered, Fig. 5 shows bifurcation plots for the same 50 different runs. All of the models do a good job matching the bifurcation point, but considerable variability is shown at long and short CLs. This result is not unexpected because data at long and short CLs were not provided for fitting. To obtain parameterizations that better reproduce these regions of the bifurcation plot, it would be necessary to include data from more CLs.

We also consider fitting all 27 parameters of the 4 V model to action potentials obtained from that model using Set 4 V-2 (human endocardial action potentials) at CLs of 500, 400, and 350 ms. Figure S2 shows the results of one run of the algorithm. The action potentials agree well at all fit CLs as well as at many unfit CLs, including 450 and 300 ms. Disagreement occurs primarily for long CLs, where the fitting underestimates APDs, and at short CLs, where the fitting generates alternans not present in the target data set. Such discrepancies above the maximum and below the minimum CLs included when performing the fitting is not unexpected; it underscores the importance of including more data if reproducing detailed alternans properties is a priority.

Over 50 runs of the algorithm, the 4 V model allows greater flexibility in choosing parameterizations that agree well with the target data. Figure S3 shows the specific values of parameters obtained. Some parameters, such as $τso1$⁠, $τw1−$⁠, and $τw+$⁠, are largely confined to fairly small portions of the ranges within their bounds, but most parameters exhibit considerable variability. Figure S4 shows bifurcation plots obtained using the parameterizations found in the same 50 runs. Although there is good agreement at the fit CLs of 500, 400, and 350 ms, there is less consistency at longer CLs, where most parameterizations underestimate APD, and at shorter CLs, where many parameterizations produce significant alternans.

Although Figs. S2-S4 only demonstrate the use of the algorithm in fitting the human endocardial action potentials obtained from the 4 V model, we also tested fitting the 4 V model to its human epicardial action potentials, which feature a spike-and-dome morphology. The algorithm is capable of finding comparably good fits that closely reproduce the action potential shape and dynamics in that setting.

To demonstrate the flexibility of the models and the algorithm, we use them to fit microelectrode recordings obtained from cardiac tissue. Figure 6 shows results from fitting endocardial action potentials at CLs of 505, 380, and 180 ms using the 3 V model (red) and the 4 V model (blue); recall that the 4 V model is able to recover a greater range of action potential shapes than the 3 V model. Note that these are representative parameterizations obtained; many parts of the algorithm have non-deterministic components that vary over algorithm runs. Both parameterizations obtained fit the action potentials well at the CLs included in the fitting and also fit many unfit CLs well. More deviation from the experimental data is observed at the lowest CLs, where the tissue and model dynamics are both quite sensitive. Both models also slightly overestimate the voltage early in the action potential.

For the 3 V model, there is moderate variability in the parameter values obtained across 50 runs of the algorithm, as shown in Fig. 7. As was observed for model recovery (see Fig. 4), the model had no strong preference for the values for $τv1−$⁠, $τv2−$⁠, and $τw+$⁠. Noteworthy differences from the model recovery case include a preference for lower values for $τd$ and $τv+$⁠; higher values for $uc$⁠, $τo$⁠, $τsi τr$⁠, and $τw+$⁠; a much lower value for $ucsi$⁠; and a much higher value for $τw−$⁠.

Over 50 runs of the algorithm, there is much more variability for the 4 V model in the specific parameter values found in the best parameterizations, as can be seen in Fig. 8. Only a few parameters, like $τv+$ and $τw+$⁠, are highly restricted to a small region within the imposed bounds. However, there are in many cases clear differences from fitting model-derived data (see Fig. S3). For example, when fitting endocardial microelectrode data, it seems it is important to fit the values of many parameters within a narrower window, including $τw∞$⁠, $θv−$⁠, $uu$⁠, $kso$⁠, and $τv1−$⁠. Other parameters clearly require different values, such as $θv$⁠, $us$⁠, $τv+$⁠, $τso1,$ and $τw1−$⁠. Only one parameter, $τso2$⁠, obviously settled into a wider range for the endocardial microelectrode fitting, indicating that fitting microelectrode data can occur for both relatively large and relatively small values of this parameter.

The bifurcation plots produced by the 50 runs of the algorithm for each model show good agreement for most CLs, as can be seen in Fig. 9. In both cases, the most disagreement occurs for long CLs not included in the fitting (especially for the 4 V model) as well as qualitative differences with regard to the presence or absence of alternans at the shortest CLs. Because extremely short and extremely long CLs were not included in the data used to generate the fitting, it is not surprising that more variability occurs.

Similar results are obtained when fitting epicardial action potentials, as discussed in detail in the supplementary material. We also performed a statistical analysis to determine parameters that were different for each model when fitting endocardial vs. epicardial data sets. For the 3 V model fitting, five parameters had different values for endocardial and epicardial data: $uv$⁠, $uc$⁠, $ucsi$⁠, $τv+$⁠, and $τw+$⁠. When the 4 V model was used, six parameters had different values: $θv−$⁠, $w∞*$⁠, $uu$⁠, $τv+$⁠, $τso1$⁠, and $τv1−$⁠. Where the models have similar structures, some of the same parameters were identified, indicating a level of consistency across the two models in reproducing differences between endocardial and epicardial action potentials. For example, $τv+$ generally has the same role in both models, and $uv$ in the 3 V model and $θv−$ in the 4 V model both identify the threshold used in locating the step in the step-function time constant $τv−(u)$⁠.

We have shown that a genetic algorithm can be used to fit all parameter values of two flexible cardiac action potentials models to match provided input data. To find a single parameterization to fit action potentials at three different CLs, the algorithm takes 30–90 s for the 3 V model and 40–120 s for the 4 V model to run, depending on processor specifics; we obtained the faster times using 16 logical cores. Different runs of the algorithm produce different parameterizations because of its random components, and the consistency of the parameter values across multiple runs varies by parameter and model. In general, the 3 V model produces more consistent parameter values than the 4 V model, but some parameters in each model cannot be well fit, most likely because they require additional data or because the algorithm is less sensitive to that parameter compared to others.

Outside of model recovery, where the parameter values used to generate the data for fitting are known, we found that it is not useful to think of a fitting as arriving at “correct” parameter values. Indeed, our results show that most parameter values are chosen within a range that is often fairly broad. Instead, it is helpful in many cases to think of parameter values co-varying. For example, Fig. 1 shows that good results can be obtained for the 4 V model fitting for a broad range of values for $τsi$ and $τw+$ as long as they increase together following the curve shown. That is, the fitness function has only a slightly lower value for parameter pairs along the curve as they fall away from the maximum. Similarly, for the 3 V model fitting shown in Fig. 1, the increase in $τd$ (corresponding to decreased excitability) can be offset by the decrease in $uc$ (excitability threshold). In the simple two-parameter setting, the space is low-dimensional enough to allow recovery of the true maximum (or very close to it). In higher-dimensional space, the fitness function has a more complicated structure and the ways in which parameters can be varied together to achieve similar fitness values also become more complex. In the 4 V model, in particular, we have found that the algorithm can find many parameterizations that are not always very similar but achieve similar fitness values.

In addition, the consistency of fit values decreases as the number of parameters fit increases. As shown in Figs. 1 and 2, when fitting a small number of parameters, the algorithm finds the correct parameterization with a high degree of accuracy, and it does so reliably. However, as more parameters are included, the parameter space becomes more complicated, and recovering the correct values from limited data becomes more challenging. In addition, we have verified that the model is more sensitive to some parameters than to others: although it can find an accurate value for a less influential parameter when only that parameter is fit, attempts to fit the same parameter along with a single more influential parameter results in an apparent inability to fit the less influential parameter.

The two models also show different strengths. The 3 V model cannot accurately reproduce most action potential shapes, but it can still find parameterizations that result in low curve error, and the parameterizations found by the model tend to be more consistent. On the other hand, the 4 V model is capable of producing more types of action potential morphologies, but at times, it cannot find the appropriate shape because curve error still may be quite low for qualitatively incorrect shapes. In particular, it has trouble producing a spike-and-dome action potential shape when the shape changes to more triangular at shorter CLs. In addition, its parameterizations are less consistent.

Although both models can achieve good results in fitting experimental data from microelectrode recordings, some problems still can occur. For example, although the use of the biphasic stimulus reduced the likelihood of problems with the action potential upstroke, we occasionally still observed some undesirable behaviors, such as an attempt to match alternans by producing a notched upstroke. The algorithm also sometimes found a best parameterization without alternans rather than with alternans, but with the APD matching the average of the long and short alternans APDs. Further work will be needed to fine-tune the algorithm to avoid such behavior. Moreover, in the future, we intend to fit data obtained from optical mapping experiments, which will require novel ways for treating the action potential upstrokes to avoid the artificially slow upstrokes that would result from directly matching the spatially averaged data represented by optical action potentials.

Our experiences in finding parameterizations to match various data sets has led us to new ways of understanding the 3 V and 4 V models and their parameters. First, we expected that not all parameters could be fit using the data provided. For example, in the 3 V model, the parameter $τv2−$ associated with the conduction velocity restitution curve slope in tissue could not be assessed from single-cell data, and the parameter $τv1−$ associated with the minimum diastolic interval could not be assessed using only data from a small number of CLs. These expectations were borne out in the results, which showed significant variability in fitting these parameters. Including more or different data in the fitting process, such as voltage clamp data or a stochastic stimulation protocol,²⁶ may be helpful.

We also found that some parameters have a very large influence and can make fitting other parameter values, to which the models are less sensitive, quite difficult. This may be the case even when setting the algorithm to fit only a single insensitive variable allows it to do so reliably and accurately. For example, if attempting to fit only $τd$ and $τo$⁠, the algorithm would obtain good fitness values for a narrow range of $τd$ but a broad range of $τo$⁠, emphasizing the model's greater sensitivity to the value of $τd$⁠. Thus, finding parameter values can serve as another type of sensitivity analysis. Such information can be useful, but the difference in sensitivity makes fine-tuning the values of relatively insensitive parameters challenging.

For many parameters, we found a fair amount of variability in the values found, even for model recovery cases where the algorithm is attempting to find parameter values to reproduce the model's own dynamics. Overall, the 3 V model exhibits less variability in parameter values found, which suggests that its parameters are closer to being identifiable in the sense that parameters have unique input into the model based on its construction and thus can be specified. In contrast, the 4 V model shows significant variability in parameter values found, which suggests that the parameters are not identifiable—that is, different combinations of parameters can lead to similar or even indistinguishable results. For models that are not identifiable, it is not possible to identify specific parameter values needed because of interdependencies among the parameters. We expect that parameter identifiability issues will be even more common for cardiac action potential models that have hundreds of parameters: it may not be possible to obtain distinct values for each parameter because of the way in which the model was constructed.

We found that the 4 V model has difficulty achieving spike-and-dome action potential shapes when matching epicardial action potentials from microelectrode recordings. However, the same model is able to reproduce action potential data generated using its own epicardial parameter set with a consistently correct shape. These two findings suggest that either the 4 V model formulation itself is not flexible enough to match epicardial microelectrode data or that the algorithm does not settle into the part of parameter space that allows this action potential morphology. We suspect that it is the latter, as we found that spike-and-dome morphologies can be produced relatively easily when fitting only a subset of 18 4 V-model parameters. In that case, the unfit parameters likely affect the parameter space such that the parameter combinations that produce the correct morphology are much denser. Identifying suitable subsets to use in fitting may improve performance in some cases but may not provide a general solution, as they require a priori knowledge about how well the model can be fit to the system being considered.

Syed et al.¹⁸ previously used a GA to fit nine human atrial model parameters to human atrial experimental and model-derived data using a population of 100 and 100 generations. Action potentials from up to four CLs were fit simultaneously. Roulette-wheel section, single-point crossover with two children, and flip-bit mutation with a 25% probability of flipping a bit in one parameter, and standard elitism (replacement of the least fit new population member with the most fit from the previous generation) were used. Fitness was assessed as the square root of the sum of either squared error or absolute differences, and an increased weight was given to the fit at −60 mV as a surrogate for APD measurement. The bounds included a range of about 1.5–4 times the minimum value; for those parameters whose minimum value was 0 and the maximum range was 0.4. Any new parameter values that fell outside the bounds were replaced with a random value within the bounds.

Kherlopian et al.²¹ implemented a GA to fit 5 conductance parameters of the Luo-Rudy I model and 9 parameters of the model of Nygren et al. to isolated guinea pig ventricular myocytes; the fit to a human atrial model showed that a reasonable fit still could be obtained in a model with a different structure simply by rescaling conductance values. The population sizes were 20 and 40 and the numbers of generations were 20 and 60 for the models of Luo-Rudy I and Nygren et al., respectively. Pairs of successive action potentials from up to three CLs were fit simultaneously. They used tournament selection (without replacement; tournament size of 2), simulated binary crossover with a rate of 90% (the remaining 10% of the time allowed parents to pass into the next generation as their own children), and two children per pair, together with a mutation probability of 10%. Fitness was assessed using the sum of squared error, and parameter values were allowed to vary by as much as 2000% from their default settings.

Bot et al.²³ largely adopted the previous implementation scheme but with a neonatal mouse ventricular model and corresponding myocyte data. They used a population size of 40, 15 generations, and bounds of up to −90 to +200% from the default parameter values. With this configuration, it took 13–15 s to fit six conductance parameters. Groenendaal et al.²⁶ further modified the same configuration to fit 9 conductance parameters of the Faber-Rudy guinea pig ventricular model³⁷ to guinea pig ventricular myocyte data. The population was extended to 500 members and the number of generations was raised to 100, with bounds set to 0.01%–299%. They studied different types of data that could improve fitting, including the use of a stochastic simulation protocol consisting of 11 randomly timed stimuli as well as a voltage-clamp protocol. They also implemented an iterative approach to reset boundaries based on the output of previous runs. Some of these modifications may prove useful for helping to constrain parameterizations for the flexible models.

Our parameter fitting method could potentially be improved by a number of modifications. First, we consider relatively minor changes that could be made to the existing algorithm. Within selection, we could adjust tournament size slightly: larger tournaments more aggressively prune lower-fitness parameterizations, while smaller tournaments preserve diversity. For crossover, we could allow each pair of parents to produce two children, which preserves more of the characteristics of the original parents, instead of a single child, which allows traits to persist longer in the overall population. We could also modify the mutation rate and size. A different population size or number of generations could be used as well: running the GA for a large number of generations not only allows more opportunity for the algorithm to converge to the best parameterization but also requires more time. Here, we selected a fairly large population because this works well with tournament selection to explore a larger area of parameter space. We aimed to keep the computational cost low, so we balanced the population size with a small number of generations to work with our selection scheme and found good convergence within ten generations in most cases. Decreasing the bounds used for parameter values facilitates finding optimal parameterizations within the narrower bounds more quickly, but has the disadvantage of eliminating any parameterizations outside those bounds. Thus, we kept the bounds large here to consider the most general case. If more information is available about expected parameter values, decreasing the bounds could improve performance. Finally, the fitness function essentially controls what the GA converges to, so that changes to its definition will affect the quality of parameterizations found. Altering the fitness function, such as to emphasize APD without considering curve error, will produce different parameterizations that may be more suitable for some applications. We note that we performed preliminary testing of most of these modifications and found a limited effect on the quality and consistency of parameterizations obtained, with the exception of the fitness function, which as noted has a strong effect.

Other larger changes in the methods used for selection, crossover, and mutation also could be implemented and potentially could improve performance. Several options for selection exist;³⁸ the canonical approach is roulette-wheel selection, in which each parent is selected randomly with a probability equal to its normalized fitness score. We compared roulette-wheel selection and found for our case that tournament selection generally converged more quickly, but that the final fitness achieved generally was comparable between the two methods. A variety of options are also available for crossover, with the most relevant being simulated binary crossover. This method does not require conversion to a fixed-point representation and is commonly found in real-coded genetic algorithms.³⁹ For mutation, the canonical method is flip-bit, but for real values, we feel Gaussian mutation is a better choice for effecting small changes in parameterizations to explore surrounding parameter space. Another mutation option is to use an adaptive mutation strategy⁴⁰ wherein mutation is a function of fitness to produce an increased rate of mutation for parameterizations with low fitness and smaller changes in parameterizations with high fitness. The rate or size of mutation also could be varied across generations.

There are also many alternative approaches to finding parameterizations outside of GAs. Other potential optimization schemes include simulated annealing⁴¹ and evolution strategies.⁴² Traditional inverse-problem approaches also could be used, and techniques like data assimilation⁴³ can be modified to include parameter estimation.⁴⁴ Furthermore, hybrid approaches that combine different types of parameter estimation recently have been shown to be useful in applications to fitting cardiac ion current formulations.⁴⁵ 

A number of limitations are associated with our work. First, we did not consider tissue-level simulations, which can affect action potential shape, duration, and other properties like alternans.^5,7,36 However, for our single-cell simulations, we used a biphasic stimulus current exported from the electrotonic current experienced during the beginning of an action potential in a corresponding one-dimensional simulation. Although we have shown that this choice can eliminate a number of potentially poor parameterizations from consideration, it is possible that this strategy may yet be insufficient for obtaining results that apply to tissue-level simulations. Several other options exist for handling such a case, including fitting the upstroke parameters separately (not within the genetic algorithm) or keeping the bounds of the excitability parameters very tight. If these approaches do not work, the remaining solution would be to use a one-dimensional cable for fitting and accept the increased computational cost (likely around two orders of magnitude).

Although we analyzed bifurcation plots for single runs to assess accuracy in the model recovery case and across multiple runs to assess consistency, we did not look at accuracy and variability (for model recovery only) in model variables not included in the fitting (the gating variables $v$⁠, $w$⁠, and, for the 4 V model, $s$⁠). Similarly, we did not compare the magnitudes and time courses of the different currents in the model. Studying differences in these quantities could provide more insights into the quality of parameterizations obtained and the behavior of the algorithm in different settings.

As discussed before, the algorithm will not always find good fits for dynamical regimes not included in fitting, such as action potentials at long and short CLs. We also use only ten beats before comparing action potentials, which may not always be long enough for the model to reach steady state because the initial conditions are not always equally close to steady states corresponding to different parameterizations, so that it may take some parameterizations longer to reach steady state than others. A slight variation in the initial conditions (primarily by reducing the initial value of $w$⁠) can facilitate reaching steady state more quickly and also can allow the algorithm to postpone conduction block to the shortest CLs. In addition, we have used only the 3 V and 4 V models within the genetic algorithm; the performance could vary for fittings to different models. We also have fit to only a few specific types of data, so our results may not generalize straightforwardly to other types of action potential data sets (such as those belonging to mammals with very short action potentials or associated with disease states). Finally, we note that many components of genetic algorithms incorporate randomness, so that results obtained for any small number of runs may not be representative.

As our next steps, we would like to consider ways to overcome several of the limitations of the current study. First, we want to identify what types of information can be used to ensure the best possible fit over the broadest possible range of parameters—for example, how many and which CLs should be included; whether including only APD error for some CLs might lead to more robust fittings, etc. In addition, we want to evaluate the best ways to set the bounds within which parameter values are chosen, as tighter bounds improve efficiency but may close off relevant areas of the parameter space. We also intend to analyze how well the algorithm predicts tissue-level behavior and whether additional information can be used to improve such predictions. Finally, we plan to find efficient ways to extend our algorithm to work with data obtained from optical mapping recordings, which, as spatial averages, pose unique challenges in terms of action potential shape and upstroke speed.

Additional information regarding the Methods and Results is available in the supplementary material, as indicated in the text.

This work was supported in part by the National Science Foundation under Grant Nos. CMMI-124235, CMMI-1341190, and CNS-1446312.