A robust multi-objective optimization framework to capture both cellular and intercellular properties in cardiac cellular model tuning: Analyzing different regions of membrane resistance profile in parameter fitting

Summary

Mathematical models of cardiac cells have been established to broaden understanding of cardiac function. In the process of developing electrophysiological models for cardiac myocytes, precise parameter tuning is a crucial step. The membrane resistance (R_m) is an essential feature obtained from cardiac myocytes. This feature reflects intercellular coupling and affects important phenomena, such as conduction velocity, and early after-depolarizations, but it is often overlooked during the phase of parameter fitting. Thus, the traditional parameter fitting that only includes action potential (AP) waveform may yield incorrect values for R_m. In this paper, a novel multi-objective parameter fitting formulation is proposed and tested that includes different regions of the R_m profile as additional objective functions for optimization. As R_m depends on the transmembrane voltage (V_m) and exhibits singularities for some specific values of V_m, analyses are conducted to carefully select the regions of interest for the proper characterization of R_m. Non-dominated sorting genetic algorithm II is utilized to solve the proposed multi-objective optimization problem. To verify the efficacy of the proposed problem formulation, case studies and comparisons are carried out using multiple models of human cardiac ventricular cells. Results demonstrate R_m is correctly reproduced by the tuned cell models after considering the curve of R_m obtained from the late phase of repolarization and R_m value calculated in the rest phase as additional objectives. However, relative deterioration of the AP fit is observed, demonstrating trade-off among the objectives. This framework can be useful for a wide range of applications, including the parameters fitting phase of the cardiac cell model development and investigation of normal and pathological scenarios in which reproducing both cellular and intercellular properties are of great importance.

Keywords:

Simulation and modeling – Optimization – Curve fitting – Mathematical models – Action potentials – Interpolation – Peak values

Introduction

Mathematical models of cardiac electrophysiology are of the utmost importance in solving a wide range of biomedical and pharmacological problems [1]. Computational models of cardiac myocytes facilitate the understanding of important properties and phenomena at both cellular and tissue levels [2]. Cardiac myocyte models are usually represented by non-linear systems of ordinary differential equations (ODEs), in which a large number of unknown parameters need to be tuned during the processes of model development and calibration [3, 4].

For tuning the parameters, several optimization methods have been proposed to overcome the high computational burden and inaccuracy of a trial and error approach [2]. Gradient and heuristic-based approaches [5–8] are the most common methods for tuning the model parameters. In accordance with the literature review, one can notice that the optimization approaches developed for tuning cardiac cell models are mainly limited to the single objective problems of fitting action potential (AP) waveforms, ionic channel currents, or a combination thereof, i.e., properties of single cardiac cells. These optimization protocols can be inaccurate in reproducing interconnected cell properties, i.e., tissue behavior [3]. However, intercellular electrical coupling properties are critical in the context of large-scale modeling [9].

To illustrate the importance of considering intercellular features, it has been demonstrated that two distinct sets of parameters can generate identical single-cell AP waveforms [10, 11], while their intercellular characteristics are different [3]. This dissimilarity stems from the incomplete characterization of current dependencies in source-sink relationships [3]. AP propagation in the tissue occurs provided that an adequate charge is delivered to local sinks by active sources [12]. Membrane resistance (R_m) is a useful indicator to characterize this intercellular communication [13, 14].

Weidmann [15] pioneered the measurement of R_m during AP waveform and described a phenomenon called all-or-nothing repolarization (ANOR) in sheep Purkinje fiber. To demonstrate all-or-none behavior, current pulses were applied after the absolute refractory period to perturb the AP trajectory. In response to these pulses, the AP either returned to the diastolic phase (repolarized) or exhibited regenerative depolarization to follow the AP waveform. This behavior, in which a perturbation of membrane voltage during the AP leads to repolarization, is known as auto-regenerative ANOR [16, 17].

To further investigate important dynamic characteristics, including R_m and auto-regenerative ANOR, a three-dimensional (3D) visualization of time, transmembrane voltage, and current (t-V_m-I_m) was introduced by Zaniboni [17–19]. In these studies, it has been demonstrated that negative R_m [20] occurs in the late phase of repolarization. The time interval during which R_m becomes negative is called the ANOR window, and its importance in the occurrence of early after-depolarizations (EADs) has been clearly shown [17]. The membrane during this time interval behaves as a source of current, which results in surges in outward currents in response to any additional hyperpolarizing current. Similarly, depolarizing perturbations in voltage result in additional depolarizing current, which may cause EADs. Moreover, on the borders of the ANOR window, high values of R_m and nonlinearity of the I_m-V_m curve appear. In the R_m profiles of similar APs obtained with different sets of parameters in the same model or R_m profiles of different models with similar cell types, it has been demonstrated that these high R_m values can occur at different voltages [17]. Furthermore, it has been shown that measuring a high value of R_m at the end of plateau is computationally unfeasible and can be affected by discontinuities of the R_m profile [21]. However, these observations have not been directly considered in the cellular fitting approach described in this paper.

The speed of AP waveform propagation is called conduction velocity (θ). Regions of cardiac tissue with slower θ are at risk of re-entrant propagation, which is associated with arrhythmia [22]. Also, the relationship between R_m during the diastolic phase (i.e., rest phase) and θ, in a human atrial model was studied in [23], which has been demonstrated that decreasing R_m results in reduced θ. This reduction results from the increase in the electrical load on neighboring cells. Moreover, the analysis of a human ventricular model suggests that reduction of [K⁺]_o results in increased R_m during the rest phase, which can be a symptom of hypokalemia [24].

The discussion above emphasizes the benefits and rationale for including R_m as one of the objective functions during parameter fitting. In [3], Kaur et al. considered R_m values together with AP as objectives for multi-objective optimization. However, Kaur et al. [3] did not account for the fact that R_m values as a function of the transmembrane voltage, V_m, increase steeply to become singular at specific values of V_m. The dependency of the R_m profile on V_m and discontinuity of this profile pose a computational challenge for the fitting procedure that must be addressed to ensure a robust procedure.

In this paper, we propose a careful selection of regions of interest for R_m during the AP waveform. These regions are selected after the analysis of different cardiac myocyte models for the human left ventricle under different physiological and pathophysiological conditions. During the different phases of the AP, three regions are identified by excluding ranges near the points which are susceptible to singularity, i.e., when R_m goes to infinity. By specifically considering R_m in these regions, the parameter fitting process is more robust. Considering the AP, R_m waveform in the repolarization phase, and R_m value in the rest phase as objectives allows us to explore the trade-offs of objectives among cellular and intercellular properties within the framework of multi-objective optimization. Moreover, to reduce the number of optimization objectives, an interpolation approach is adapted to fit a curve to R_m points, and fitness of this curve instead of individual R_m points is used as one of the objective functions. This enhances the computational tractability of optimization and facilitates visualization of the Pareto front.

Materials and methods

Mathematical models and simulation specifications

Two mathematical models of the human ventricular epicardial cell are used in this study for the optimization: the Ten Tusscher et al. (TNNP) model [25], and the Iyer-Mazhari-Winslow (IMW) model [26]. All simulations are carried out in Matlab R2018a [27], using codes of the models available at www.cellml.org. The detailed description of ODEs associated with these human ventricular cell models can be found in [25, 26]. The built-in ode15s solver in Matlab for stiff systems of ODEs is utilized. Also, the Parallel Computing Toolbox is used in the simulations to reduce the computational time.

The AP is produced by injecting a 2 ms current pulse with 20 pA/pF amplitude [10]. To reach the steady-state condition, each model is allowed to run for 20 beats at a 1 s interval (1 Hz), and the 20^th AP is selected for analysis. Values of R_m are normalized to a cell capacitance of 180 pF [17, 18].

The TNNP and different configurations of IMW models are considered as the base and reference (experimental data) models, respectively. In other words, key parameters of the base model (TNNP) are adjusted to attempt to reproduce data generated by the reference model (IMW). The importance of selected parameters is thoroughly discussed in [10]. In addition to the mentioned models, another state-of-the-art ventricular model [28] is utilized to further examine the R_m profile for the proper selection of regions of R_m far from singular points.

R_m measurement

The protocol used to measure R_m during the AP waveform at a particular point (voltage) during repolarization is illustrated in Fig 1. For each specific V_m value of interest, the ODEs of the model are solved in two separate simulations. Both simulations start by simulating an AP until reaching the point (voltage) of interest during repolarization. At this point, the simulation is changed to voltage-clamp mode, and the voltage is kept constant at 10 mV higher (V_m+10) or lower (V_m-10) than the point of interest, respectively. The resulting membrane currents (I_m+10 and I_m-10) are recorded after 5 ms of voltage-clamp in each simulation. The ratio of difference in voltage (ΔV_m) to change in current (ΔI_m) is used to estimate R_m [3, 14].

To avoid ambiguities where the same voltage value occurs more than once during repolarization, the AP waveform is divided into three parts based on global/local minima and maxima (red marks in Fig 1). In each part, the voltage clamp protocol described above is applied. Fig 1 represents the R_m profile (brown line) calculated from the AP waveform of the TNNP model as an example. In the R_m profile, the ANOR window (grey rectangle) [17] is bounded by the sharp increase of positive and negative R_m values. Inside the window, R_m values are negative.

It is worth noting that the R_m profile is constructed as a function of voltage, rather than specific points in time. This allows a more direct comparison between models even when there may be a substantial discrepancy between models in the timing of different states of the AP (i.e., substantial differences in AP shape and duration) [3].

Discontinuity and voltage dependence of the R_m profile

As illustrated in Fig 1, R_m is positive during the plateau phase of the AP, becomes negative during the ANOR window, and returns to positive values during late repolarization. The transition from positive to negative values (and the reverse) is associated with a discontinuity due to singular R_m points, i.e., the magnitude of R_m increases, rapidly approaching infinity. Fig 2 provides further insight into the rapidly changing R_m near this discontinuity. In this figure, the R_m profile of the TNNP model is calculated by applying voltage clamp using a range of sampling resolutions (1 mV, 0.5 mV, 0.1 mV). Fig 2 confirms that decreasing the step size to record high resolution of R_m profiles results in increasing peak values for R_m (red oval). Therefore, depending on the R_m step size, different peak values of R_m will be calculated (approaching infinity as the voltage step size decreases). Hence, attempts to match the R_m profile for parameter tuning near these transition points cannot be expected to yield robust results. However, as can be perceived from the parts of R_m profile highlighted by green ovals in Fig 2, this is not a concern other than in the immediate vicinity of the discontinuities.

The challenge described above is further compounded by the fact that R_m points, which are susceptible to be singular, appear at different values of V_m for different cell models, even when the cell models have similar AP (Fig 3(A)). This can be observed in Fig 3(B), which shows that the peak values of the R_m profile of the fitted base model (fitting result) are not aligned with the peak values of the reference model (IMW).

Based on the above discussion, it is important to carefully select the best regions of R_m to use for parameter fitting. This is performed by examining different cell models in different states, as described next.

Defining allowed and disallowed regions for capturing R_m in optimization

To establish specific disallowed regions for our models of interest, for each model configuration of interest, we first find the singular points (see example in Fig 4). We then exclude voltage ranges containing these points in any model of interest (“Disallowed Regions” in Fig 4) from consideration and focus the fitting of R_m on the voltage ranges with no singularities (“Allowed Regions” in Fig 4). In the proposed approach, objective functions related to R_m are selected from allowed regions only, eliminating the discontinuity problem discussed in the previous subsection.

To define the allowed and disallowed regions for this paper, three human cardiac ventricular cell models (TNNP, different configurations of the IMW, and [28]) are used. The extent of the disallowed regions is chosen as 5 mV outside the “outermost” singularity (black pentagrams in Fig 4) of any of the models of interest. The singularity in repolarization phase for one of the considered ventricular models [28] (not shown in Fig 4) happens around -70 mV; as a result, the extension of Disallowed Region 2 is similar to what is seen in Fig 4. The determination of the disallowed regions is specific to the models and AP configurations of interest in a particular fitting problem and must thus be carried out for the defined fitting scenario of interest.

Problem definition of the multi-objective optimization

The general formulation of the multi-optimization problem (MOP) is as follows:

F(x) reflects cellular and intercellular parameters. Amid cellular properties, minimizing the root-mean-squared-error (RMSE) of AP is considered as the first objective. The RMSE of AP (RMSE_AP) can be calculated as follows:

Fitting selected points of the R_m profile in the allowed regions, discussed in the previous subsection, could be taken into consideration as a key intercellular property. In this study, errors in fitting several points within Allowed Region 2 are considered as potential objectives of the optimization problem. However, considering multiple resistance points increases the dimension of the optimization problem, thereby increasing the computational complexity. Hence, we instead settle on the curve of R_m profile obtained in Allowed Region 2 as the second objective. To calculate this curve, first, the voltage clamp protocol described in R_m measurement subsection is applied to the AP waveform within Allowed Region 2 for each 5 mV. Then, the calculated discrete points of R_m are connected by interpolation to produce a more densely sampled R_c curve. The RMSE of R_c (RMSE_Rc) as a second objective is then obtained as follows:

The value of R_m in the diastolic (resting) phase (R_d) in Allowed Region 3 is another intercellular property derived from the R_m profile. The absolute error (AE) of R_d (AE_Rd) is utilized as a third objective:

Based on the above description, three objective functions are introduced, in which M varies regarding the scenarios defined in the description of optimization scenarios subsection.

The x_i (i = 1,…,16) are G_Na, G_CaL, G_bNa, G_bCa, G_Kr, G_to, G_K1, G_Ks, G_pk, P_NaK, k_NaCa, G_pCa, V_maxup, c_rel, a_rel, and V_leak for the results presented in this paper. These parameters are maximum conductances and maximum rates for transporters, pumps, and intercellular Ca²⁺ handling mechanisms in the model. The ub_i and lb_i are defined within the physiological ranges and obtained by trial and error search [6]. More explanations of parameters and boundaries are provided in Table A1 in the S1 Appendix. Our formulation does not require equality and inequality constraints described in Eq 1.

The solution to the multi-objective problem

The outcome of the optimization is a set of optimal solutions, in which there is a trade-off among different conflicting objectives. In other words, there are contradictions among objectives, and some of them can only be optimized while others are not at their optimum value. The basic terminologies regarding MOP solutions are represented as follows and are further explored in Fig 5.

Definition 1- Dominant solution: A solution set x₂ dominates x^*, if f_i (x₂) ≤ f_i (x^*), i = 1, 2,…, M (where M is the number of objectives) and among the different inequalities, at least one of them is strict.

Definition 2- Pareto front: A solution that is not dominated by other solutions is a Pareto-optimal solution. The set of all Pareto solutions constitute the Pareto front.

In Fig 5, x₁, x₂, and x₃ are Pareto optimal solutions; and therefore, they are among the solution sets which held the trade-off among objectives. As can be seen in this figure, Pareto optimal solutions from the Pareto front (black curve), which consisting of the solutions of the Pareto front set, are not dominated by other solutions. In the case that the Pareto front is obtained precisely, no feasible solution below the Pareto front can be attained. Therefore, the solution point, marked with yellow, is infeasible. The green point called utopia is the minimum solution for all the objectives; however, reaching this solution is infeasible. It is worth noting that the multi-objective optimization aims at attaining solution sets, forming the Pareto front approximation with minimum deviation from the true Pareto front.

Multi-objective evolutionary algorithms background

The most general and simplest approach to solve the MOPs is the weighted-sum approach, in which the user defines the values of weights. The choice of weight values highly depends on the problems and the level of priority of each objective. In this method, the MOP is transferred to a single objective by defining a new objective function, which is the summation of each normalized weighted objective [29]. Being dependent on the weight vector, selected by trial and error in a time-consuming procedure, is amid the major drawbacks of this approach. This problem can be addressed by multi-objective evolutionary algorithms (MOEAs) in which the simultaneous search of algorithm results in finding the trade-off among objectives. A wide range of MOEAs has been proposed in the last decade to solve complex optimization problems heuristically [30, 31] in various fields of study, including biological engineering [32–34].

NSGA-II

Amongst MOEAs, NSGA-II (non-dominated sorting genetic algorithm II) [35] is a powerful optimization tool. Some properties of this algorithm for utilizing in our problem are as follows:

1) Utilizing the elitism principle maintains the best solutions for the next generation, 2) crowding distance technique provides a diversity of the solutions, and 3) it has been demonstrated to be a successful algorithm in many engineering problems.

The definitions required to understand the procedure of NSGA-II are given as follows:

Definition 3- Crossover: The crossover is an operator to mix the genetic information of the parents and produce the Offspring. In NSGA-II, simulated binary crossover (SBX) is used, which creates pairs of offspring from parents based on the crossover index.

Definition 4- Mutation: The genetic diversity of the population from one generation to another is controlled by the operator known as mutation. In NSGA-II, the polynomial mutation is utilized, which controls the probability distribution of solutions using an external parameter. This parameter adjusts the diversity of solution sets.

Definition 5- Crowding distance: Crowding distance is determined by the cube enclosing each solution and measures the diversity of a solution from its neighbor. The solution with the largest crowding distance wins the tournament selection and is kept for the next iteration.

A brief description of NSGA-II procedure is as follows:

Step 1- Initialization: at the initial iteration, t = 0, the parent population of size N (P_t) is randomly generated by assigning random values to the parameters within the upper and lower bounds as described in the problem definition of the multi-objective optimization subsection. Crossover and mutation operations are utilized to generate the offspring population (Q_t) of the first iteration.

Step 2- Pareto fronts construction: to extract the Pareto fronts (F_i, i = 1, …), firstly, the combined population (R_t) consisting of P_t and Q_t is generated at each iteration. Then, non-dominated sorting is applied to each solution of the population to categorize the R_t into fronts. Individuals who are not dominated by others are chosen for the 1^st front. Such individuals that are dominated by individuals in the 1^st front, while they are not dominated by any other individuals, are classified as front 2. This procedure is repeated recursively to form all possible Pareto fronts.

Step 3- Elitist set construction: the individuals participate in the competition to be selected. The criteria to choose the elitist sets of parameters (problem definition of the multi-objective optimization subsection) are based on the rank of domination and crowding distance. For individuals in different fronts, those in the lower front win the tournament (non-dominated rank). If individuals have the same rank (same front), the individuals with larger crowding distance are selected. The elitist set is used as the parent population of the t+1^th iteration (P_t+1).

Step 4- Satisfying the stopping criterion: mutation, crossover, and steps 2–3 are repeated until the stopping criterion is met. Here, a maximum number of iterations is set as the stopping criterion.

Finally, the 1^st Pareto front of R_t obtained from the last iteration is attained as the solution sets.

As population size and the number of evaluations greatly impact the execution time of the optimization, by performing a sensitivity analysis using a grid search and comparing the Pareto sets of various runs, 100, 10000, 20, and 20 are opted for population size, the number of evaluations, the distribution index of simulated binary crossover, and distribution index of polynomial mutation, respectively. Therefore, the results of runs pertain to this set of hyperparameters are reported in this paper.

Selection of optimal solution

Once the Pareto front for a multi-objective optimization problem has been found, the selection of a single optimal solution as a preferred solution is essential in the decision-making process. In this work, the following equation is defined to choose a desirable solution:

Results

Description of considered AP configurations

In this paper, three scenarios consisting of different combinations of objectives (see next subsection) are utilized to examine the efficacy of the proposed optimization problem in fitting the characteristics of the base model to the reference one. In each scenario, characteristics of five distinct configurations of the IMW model (Fig 6) are considered as references to verify the robustness of the approach. In Configuration 1, the baseline IMW model is used as a reference. Configurations 2 and 3 explore parameter sets that produce shortened and prolonged action potential duration (APD) of the IMW model, respectively. These configurations are defined based on the crucial role of potassium (K⁺) channels in arrhythmia disease. Increased K⁺ current amplitude results in APD shortening and can induce re-entry, while decreased K⁺ current amplitude prolongs APD. EAD and a form of abnormal heart rhythm called torsades de pointes can occur due to excessive prolongation of the APD [37]. To simulate these conditions, the availability of delayed rectifier K⁺ currents (I_kr and I_ks) are changed by manipulating the maximum conductances associated with these currents. Configurations 4 and 5 are devised to represent conditions of hyperkalemia and hypokalemia, respectively. In hyperkalemia [K⁺]_o is increased to above 6 mM, and in hypokalemia [K⁺]_o is decreased to below 3.5 mM [24]. Thus, in Configurations 4 and 5, [K⁺]_o are set to 6.5 mM and 3 mM, respectively.

As described in the problem definition of the multi-objective optimization, the RMSE and AE, well-defined evaluation criteria in the literature (e.g., in [5]), are used to assess the accuracy of the fitting.

Description of optimization scenarios

Three distinct optimization scenarios are considered in this paper to assess the contributions of different fitting objective functions to the outcome of the optimization.

Scenario 1 involves single-objective optimization considering only AP fitness. The optimization is defined and solved using GA. In this scenario, RMSE_AP obtained by Eq 2 is the only objective function.

Scenario 2 is a bi-objective optimization problem with two objective functions, namely RMSE_AP (Eq 2) and RMSE_Rc (Eq 3). The optimization is solved by NSGA-II.

In Scenario 3, minimizing the AE_Rd (Eq 4) is considered as an additional objective in addition to reducing the RMSE_AP and RMSE_Rc. The multi-objective optimization is solved by NSGA-II.

The relative percentage of improvement or deterioration in fit for the different scenarios is calculated as follows:

For the purpose of comparing the accuracy of various scenarios in fitting RMSE_Rc and AE_Rd, a posteriori calculations of RMSE_Rc and AE_Rd are conducted for Scenario 1, while a posteriori calculation of AE_Rd is performed for Scenario 2.

Numerical results

The average (Ave.) and standard deviation (Std.) of RMSE_AP (mV), RMSE_Rc (GΩ), and AE_Rd (GΩ) in 10 trials are calculated for all scenarios and configurations, and the results are presented in Table 1. In Scenario 1, the range of Ave. and Std. of RMSE_AP are changed within [1.05, 1.88] mV and [0.15, 0.28] mV, respectively. Although, the AP fit in this scenario is relatively good (Fig 6), RMSE_Rc is quite high. As an example, this poor-fitting R_c in Scenario 1 is demonstrated for one of the configurations in Fig 7. This figure illustrates the difference between R_c of the optimization result and the reference for one of the trials of Scenario 1, Configuration 3 (RMSE_AP = 1.53 mV, RMSE_Rc = 0.35 GΩ). This observation demonstrates that considering the AP only as an objective may not necessarily result in a good fit for intercellular properties (i.e., R_c) to those of the reference (model or experimental data). As can be clearly seen from Table 1, in Scenario 2, the Ave. of RMSE_Rc in Configurations 1–5 are considerably improved by 89.89%, 95.05%, 88.34%, 95.52%, and 96.14%, respectively in comparison to corresponding configurations in Scenario 1. Fig 7 confirms this significant improvement in RMSE_Rc compared to Scenario 1. Moreover, Std. of RMSE_Rc in Scenario 1 are notably larger than the corresponding values in Scenario 2 (Table 1). This observation shows that fitting AP only (Scenario 1) does not guarantee the consistent fitting of R_m, and generally results in high variability in R_m fitting. Scenario 2 results in a robust fitting of R_m, as evidenced by a significant improvement in the Ave. and Std. of RMSE_Rc compared to Scenario 1. However, this comes at the expense of relatively worse results for the AP fit. Specifically, the Ave of RMSE_AP increases by 25.93%, 33.07%, 22.76%, 20.49%, and 30.40% in Configurations 1–5, respectively, compared to the corresponding values for Scenario 1.

As shown in Table 1, in Scenario 3, minimization of the third objective (AE_Rd) is realized at the expense of increasing Ave. of RMSE_AP by 60.34%, 69.55%, 52.87%, 40.82%, and 13.92% in Configurations 1–5, respectively, compared to Scenario 2. In this scenario, the Ave. of RMSE_Rc is better than in Scenario 1, and the Std. of RMSE_Rc is slightly worse than in Scenario 2. Importantly, depending on the specific configuration, Scenario 3 can result in up to 98.58% improvement (Configuration 5) in Ave. of AE_Rd compared to two other scenarios. Compared to Scenario 1 and 2, Scenario 3 yields slightly lower Std. of AE_Rd for all configurations.

Three trials of Scenario 1 for Configuration 1 regarding AP with similar RMSE are shown in Table 2 to compare their RMSE_Rc. As can be discerned from Table 2, the fitting results in terms of RMSE_AP are comparable for all three trials; however, there is significant variability in RMSE_Rc. Hence, this table further demonstrates that different sets of parameters can yield similar cellular property (i.e., AP), while their R_m values—indicators of intercellular property—are quite different.

To further examine the performance of the proposed approach, the Pareto fronts of all configurations for Scenario 2 are shown in Fig 8. Fig 8A–8E represents examples of the convex Pareto fronts attained from solving Scenario 2 for different configurations. It is clear from these figures that there is a trade-off between minimizing the RMSE_AP and RMSE_Rc. The selection of the optimal solution from the Pareto front is highly dependent on the application and whether AP or R_c improvement is considered more important. In this general study, relying on existing literature on the applications of multi-objective optimization, we select one preferable solution as described in the selection of optimal solution subsection. It can be seen from Fig 8 that, in all configurations, the selected points yield a lower error in comparison to Ave. of RMSE_Rc in Scenario 1. Hence, these results also demonstrate that RMSE_Rc can be improved only by sacrificing AP fitness.

Fig 9 shows one trial of Configuration 2 in Scenario 3, illustrating the trade-off among objectives. The Pareto front obtained for this configuration is shown in Fig 9(A). In this figure, the selected solution set is indicated by an arrow (RMSE_AP = 2.45 mV, RMSE_Rc = 0.018 GΩ, and AE_Rd = 0.0005 GΩ) and the color scale is associated with the third objective (AE_Rd). To further elaborate on trade-offs among objectives, the heatmap of Pareto front is demonstrated in Fig 9(B). Solutions are ordered based on the low value to the high value of the first objective (RMSE_AP). The color of each cell shows the normalized error for the corresponding objective. The confliction of objectives can be vividly interpreted in this figure.

Discussion

Comparison with previous work

In this paper, a novel problem formulation for fitting human ventricular cellular models is proposed. The work presented in [3] was an important starting point and inspiration for our work. In [3], 9 parameters were optimized, and a small number of arbitrarily selected R_m points considered in a model-to-model fit. To obtain more accuracy, we have explored a wider range of search space by considering 16 parameters, defined in the problem definition of the multi-objective optimization, for fitting similar models to those used in [3]. The study in [3] found that AP only fitting did not produce accurate AP fits, and that also considering R_m resulted in improvement of the AP fit. However, with more investigation in this study, it is observed that acceptable AP fitting can be achieved using the error in the AP waveform as the only objective function with higher search space dimensions (Fig 6). Here, the proposed approach is also evaluated using four other configurations of the IMW model, including shortened APD, prolonged APD, hyperkalemia, and hypokalemia conditions.

It is important to note that in Scenario 2, an arbitrary selection of R_m values as fitness functions may be susceptible to choosing the discontinuous R_m points. As an improvement to the cellular model tuning proposed in [3], we define allowed regions within which R_m values can be used as objective functions and avoid regions near the singularities of R_m. Thus, the R_m points are selected from regions that are less vulnerable to producing errors in model fitting. Moreover, in this work, the number of objectives is reduced by curve fitting instead of fitting the distinct points of R_m. This reduction simplifies the optimization and yields enhanced decision-making ability on the final solution.

Trade-offs among objectives

In this work, we have fitted up to three different objective functions, RMSE_AP (mV), RMSE_Rc (GΩ), and AE_Rd (GΩ). The RMS error in matching the resultant curve is considered as a second fitting objective function. The main reason to select the R_m curve in Allowed Region 2 (RMSE_Rc) as one of the objective functions is that negative values of R_m and ANOR phenomenon emerge in this voltage range [17–19]. The convex Pareto-front in the bi-objective optimization problem solved by NSGA-II is demonstrated (Fig 8) to select a preferred solution based on the application. Constraining the fitting problem by increasing the number of objective functions (i.e., RMSE_Rc) results in a substantial reduction of average Std. of RMSE_Rc fitting (Table 1). This reduction demonstrates the repeatability of the algorithm across trials. The importance of R_d in determining the conduction velocity θ–one of the intercellular properties–and its relationship with [K⁺]_o is discussed in [23, 24]. This is a justification for the importance of R_d and consideration of AE_Rd as an objective in the proposed problem formulation. Our observations reveal that there tends to be a trade-off between this property and the two other features considered, i.e., RMSE_AP and RMSE_Rc. Fig 9 shows that improvements in terms of minimizing R_d are feasible by sacrificing the accuracy in fitting AP and/or R_c. However, the Ave. and Std. of RMSE_Rc in this scenario is still improved in comparison with Scenario 1.

Further investigation is required to clarify the best decision-making scheme for selecting the final solution from the Pareto front. The careful selection of additional voltage-dependent parameters and their consideration in the optimization may result in further improvement in the fitting results as it may provide an opportunity to better align the voltage-dependence of the R_m profile with the reference data. However, additional experiments that are beyond the scope of this study would be required to determine whether this would indeed improve overall fitting results.

Future work

A key limitation of the current work is the requirement for the selection of allowed and disallowed regions, which depend on the specific models and configurations used. As a future development, it would be desirable to develop a generalized algorithm to eliminate the dependency on the designers’ somewhat subjective decisions in this regard.

The model fitting can involve some degree of uncertainty due to the variability of observations. Populations of models (POMs) are used to model the inherent variability of experimental data [38]. As a further study, the proposed framework for R_m could be extended to POM calibration. Considering R_m features might reduce the degree of uncertainty in POM as it can provide intercellular information.

Conclusion

In this study, the benefits of including R_m values in cardiac cellular fitting are clearly demonstrated. We have also identified important technical challenges related to the occurrence of discontinuities in R_m and the voltage dependence of the R_m profile. These need to be addressed in the formulation of the optimization problem, or significant errors may result. This work presents and evaluates a new problem formulation for tuning specified parameters of cell models. Importantly, we provide information about the trade-offs among fitting objectives and demonstrate that it is possible to achieve a good balance among key cellular and intercellular properties. In conclusion, this paper contributes toward making mathematical models of human ventricular cells (and cardiac cells in general) more reliable by developing a robust optimization formulation for cardiac cellular fitting.

Supporting information

S1 Appendix [docx]
Overview of cardiac cellular model parameters.

S1 Text [docx]
Description of data.

S1 File [zip]
Data.