Gradient-based parameter optimization method to determine membrane ionic current composition in human induced pluripotent stem cell-derived cardiomyocytes

Premature cardiac myocytes derived from human induced pluripotent stem cells (hiPSC-CMs) show heterogeneous action potentials (APs), probably due to different expression patterns of membrane ionic currents. We developed a method for determining expression patterns of functional channels in terms of whole-cell ionic conductance (Gx) using individual spontaneous AP configurations. It has been suggested that apparently identical AP configurations can be obtained using different sets of ionic currents in mathematical models of cardiac membrane excitation. If so, the inverse problem of Gx estimation might not be solved. We computationally tested the feasibility of the gradient-based optimization method. For a realistic examination, conventional 'cell-specific models' were prepared by superimposing the model output of AP on each experimental AP recorded by conventional manual adjustment of Gxs of the baseline model. Gxs of 4–6 major ionic currents of the 'cell-specific models' were randomized within a range of ± 5–15% and used as an initial parameter set for the gradient-based automatic Gxs recovery by decreasing the mean square error (MSE) between the target and model output. Plotting all data points of the MSE–Gx relationship during optimization revealed progressive convergence of the randomized population of Gxs to the original value of the cell-specific model with decreasing MSE. The absence of any other local minimum in the global search space was confirmed by mapping the MSE by randomizing Gxs over a range of 0.1–10 times the control. No additional local minimum MSE was obvious in the whole parameter space, in addition to the global minimum of MSE at the default model parameter.

Step size to move NP x Subscript to represent membrane current, such as I Na , I CaL , I K1, I ha , I Kr , I Kur , I Ks , and I bNSC Over the past half-century, the biophysical characteristics of ion-transporting molecules (channels and ion exchangers) have been extensively analyzed. Biophysical models of each functional component have largely been detailed [1][2][3][4] , including human induced pluripotent stem cells (hiPSC-CMs) [5][6][7] . In addition, various composite cell models, including membrane excitation, cell contraction, and intracellular ionic composition homeostasis, have been developed by integrating mathematical models at the molecular level into cardiac cell models [8][9][10][11] . These models have been useful for visualizing individual currents underlying the action potential (AP) configuration under various experimental conditions in mature cardiac myocytes. However, the utility of these mathematical cell models is limited because of the lack of extensive validation of the model output accuracy. This is a drawback of the subjective manual fitting method used in almost all published mathematical cardiac cell models. A new challenge of mechanistic models of cardiac membrane excitation might be an examination in a very different paradigm to assess if the many, but continuous, variety of cardiac AP configurations, such as those recorded in hiPSC-CMs, can be reconstructed by applying the automatic parameter optimization method to the hiPSC-CM version of human cardiac cell models. We do not intend to propose a new hiPSC-CM model. The automatic parameter optimization technique objectively determines parameters in a wide range of biological models, including cardiac electrophysiology [12][13][14][15] , systems pharmacology [16][17][18][19][20] , and other models. Because of this utility, many improvements in information technology have been realized 21,22 . However, in electrophysiology, different combinations of model parameters may produce very similar APs 13,[23][24][25] . The determination of current density at high fidelity and accuracy likely requires additional improvements to the optimization method in the cardiac cell model because of the complex interactions among ionic currents underlying membrane excitation 23,26 .
The final goal of our study is to develop an objective and accurate method for determining the current profile (i.e., the expression level of functional ionic currents) underlying individual AP configurations. As a case study, we chose a large variety of AP configurations in hiPSC-CMs, which are difficult to classify into the conventional nodal, atrial, or ventricular types. The molecular bases of the ion channels expressed in hiPSC-CMs correspond to those in adult cardiac myocytes in the GSE154580 Gene Expression Omnibus Accession viewer. Electrophysiological findings suggest that the gating of ionic currents is quite similar to that observed in mature myocytes 27 . Thus, we modified the ion channel gating kinetics of the human ventricular cell (hVC) model 11 according to the prior experimental measurements 27 for a hiPSC-CM type baseline model of the parameter optimization (PO) method. For simplicity, we assumed that the opening/closing kinetics of ion channels expressed by the same human genome remains the same among hiPSC-CMs. We also assumed that the heterogeneity of the electrical activities of hiPSC-CMs might be determined by the variable expression levels of ion channels in the cell membrane. We computationally examined the feasibility of one of the basic gradient-based optimization methods, the pattern search (PS) algorithm 21,22,28 , in a model of cardiac AP generation. We prepared a given AP configuration using each ' cell-specific model' prepared by the conventional manual fitting of the hVC model to the respective experimental recordings. To assess the accuracy of the PS method for parameter optimization, the AP waveform generated by the cell-specific model was used as a target of the optimization. The initial set of parameters for the optimization was then prepared by uniform randomization centered around the default values of the model. The PS algorithm should return the original parameter values by decreasing the mean squared error (MSE) function between the modified model output and target AP waveforms. The accuracy of the optimization was determined by recovering the original values of each ionic current amplitude as the MSE progressively decreased toward zero.

Materials and Methods
Baseline model of hiPSC-CM membrane excitation. The baseline model of hiPSC-CMs was essentially the same as the hVC model, which has been detailed 10,11 and which shares many comparable characteristics with other published human models 8,9 . The hVC model consists of a cell membrane with a number of ionic channel species and a few ion transporters, the sarcoplasmic reticulum equipped with the Ca 2+ pump (SERCA), and the refined Ca 2+ releasing units coupled with the L-type Ca 2+ channels on the cell membrane at the nanoscale dyadic space 4,29 , contractile fibers, and three cytosolic Ca 2+ diffusion spaces containing several Ca 2+ -binding proteins (Fig. S1). All model equations and abbreviations are described in the Supplemental Materials.
The kinetics of the ionic currents in the baseline model were readjusted according to new experimental measurements if available in hiPSC-CMs 27 (Fig. S2). In the present study, the net membrane current (I tot_cell ) was calculated as the sum of nine ion channel currents and two ion transporters (I NaK and I NCX ) (Eq. 1).
The membrane excitation of the model is generated by charging and discharging the membrane capacitance (C m ) using the net ionic current (I tot_cell ) across the cell membrane (Eq. 1). The driving force for the ionic current is given by the potential difference between V m and the equilibrium potential (E x ) (Eq. 3). The net conductance of the channel is changed by the dynamic changes in the open probability (pO) of the channel, which is (1) I tot_cell = I Na + I CaL + I ha + I K1 + I Kr + I Ks + I Kur + I Kto + I bNSC + I NaK + I NCX The MSE was stabilized by obtaining a quasi-stable rhythm of spontaneous APs through continuous numerical integration of the model. Typically, 30-100 spontaneous cycles were calculated for a new set of sf x s. The MSE was calculated within the time window. The width of the time window was adjusted according to the AP phase of interest. where N is the number of digitized V m points with a time interval of 0.1 ms.
In typical parameter optimization, V m,a is generated by modifying the baseline model for comparison with the experimental record (V m,t = V m,rec ). However, to evaluate the identifiability of the parameter optimization, a simple approach was adopted in the present study. We used the manually adjusted ' cell-specific' model for the target (V m,t ), which was nearly identical to V m,rec . More importantly, the ' cell-specific' V m is completely free from extra fluctuations (noise), which were observed in almost all AP recordings in hiPSC-CMs. In the optimization process, the initial value of each optimization parameter was prepared by randomizing the sf x s of the cell-specific model by ± 5-15% at the beginning of each run of PS (V m,orp ) in Eq. 8, and several hundred PS runs were repeated. Thus, the error function is This optimization method was termed the ' orp test' in the present study. The advantage of using a manually adjusted cell model for the optimization target is that the accuracy of parameter optimization is proved by recovering all sf x = 1 (Eq. 6) independent from the randomized initial parameter set. The same approach was used in a previous study 23 to evaluate the accuracy of parameter optimization by applying the genetic algorithm to the TNNP model of the human ventricular cell 33 .
Optimization using the randomized initial model parameters was repeated for more than 200 runs. Thus, the orp test might be classified in a 'multi-run optimization' . The distribution of the sf x data points obtained during all test runs was plotted in a single sf x -MSE coordinate to examine the convergence of individual sf x s with the progress of the orp test.
Pattern Search method for parameter optimization. For a system showing a relatively simple gradient of MSE along the parameter axis, gradient-based optimization methods are generally more efficient than stochastic methods for this type of objective function. We used the Pattern Search (PS) algorithm, a basic gradient-based optimization method. The computer program code for the PS 34 is simple (see Supplemental Materials) and does not require derivatives of the objective function. We implemented the code in a homemade program for data analysis (in VB) to improve the method for better resolution and to save computation time.
The primary PS method uses the base and new points 28 . Briefly, sf x is coded with symbols BP x and NP x in the computer program, representing a base point (BP x ) and a new search point (NP x ), respectively. The MSE is calculated for each movement of NP x by adding or subtracting a given step size (stp) to the BP x , and the search direction is determined by the smaller MSE. Then, the entire mathematical model is numerically integrated (Eq. [2][3][4][5] using NP x to reconstruct the time course of AP (V m,a ). This adjustment is performed sequentially for each of the four to six selected currents in a single optimization cycle. The cycle is repeated until no improvement www.nature.com/scientificreports/ in the MSE was achieved by a new set of NP x s. The BP x set is then renewed by the new set of NP x for the subsequent series of optimizations. Simultaneously, stp is reduced by a given reduction factor (redFct of 1/4). The individual PS run is continued until the new stp became smaller than the critical stp (crtstp), which was set to 2-10 × 10 -5 in the present study.

Selection of ionic currents for the optimization.
When we obtain a new experimental record of AP, we do not start the analysis with an automatic optimization of G x . Rather, we first adjust the baseline model by conducting conventional manual fitting. The nine ionic currents in Eq. 1 in the baseline model are adjusted incrementally to superimpose the simulated AP on the experimental one. During this step, it is important to pay attention to the influence of each sf x adjustment on the simulated AP configuration on the computer display. Thereby, one may find several key current components that should be used for the automatic parameter optimization. Usually, currents showing a relatively large magnitude of G x were selected for automatic optimization according to Eq. 1, while those that scarcely modified the simulated AP were left as default values in the baseline model.
Principal component analysis of cell-specific models. When the orp test was performed with p elements, it was possible to record the final point BP, where the MSE was improved in the p-dimensional space. Suppose that we represent the matrix when n data points are acquired as an n × p matrix X. In that case, we obtain a vector space based on the unit vector that maximizes the variance (first principal component: PC1) and the p-dimensional unit vector orthogonal to it (loadings vector w (k) = w 1 , w 2 · · · , w p ) . It is possible to convert each row x (i) of the data matrix X into a vector of principal component scores t (i) . The transformation is defined as To maximize the variance, the first weight vector w (1) corresponding to the first principal component must satisfy: The kth component can be determined by subtracting the first (k-1)-th principal components from X: The weight vector is then given as a vector such that the variance of the principal component scores is maximized for the new data matrix: Membrane excitation and its cooperativity with intracellular ionic dynamics. When any of the G x s is modified, the intracellular ion concentrations ([ion] i ) change, although the variation is largely compensated for with time in intact cells by modifying the activities of both the three Na + / two K + pump (NaK) and three Na + / one Ca 2+ exchange (NCX). In the present study, we imitated the long-term physiological homeostasis of [ion] i by introducing empirical Eq. 13 and 14. These equations induced 'negative feedback' to the capacity (maxI NaK and maxI NCX ) of these ion transporters. Each correcting factor (crf x ) was continuously scaled to modify the limiting activity of the transporters to maintain the [Na + ] i or the total amount of Ca within the cell (Ca tot ) equal to their pre-set level (std Nai , std Catot ) with an appropriate delay (coefficients 0.3 and 0.008 in Eq. 13 and 14, respectively).
For the control of [Na + ] i : For the control of Ca tot : Ca tot is given by [Ca] i included in the cytosolic three Ca-spaces jnc, iz, and blk, and in the sarcoplasmic reticulum SR up and SR rl in the free or bound forms, respectively.
where vol is the volume of the cellular Ca compartment (see more details 11 ).
Preparation of dissociated hiPSC-CMs and recording of spontaneous APs. The 201B7 and 253G1 hiPSC lines generated from healthy individuals were used in this study 35,36 . The differentiation of hiPSCs www.nature.com/scientificreports/ into cardiomyocytes was promoted using an embryoid body differentiation system 37 . The hiPSCs were incubated at 37 °C in 5% CO 2 , 5% O 2 , and 90% N 2 for the first 12 days to promote differentiation. hiPSCs aggregated to form embryoid bodies and were cultured in suspension for 20 days. On day 20 of culture, embryoid bodies were treated with collagenase B (Roche, Basel, Switzerland) and trypsin-EDTA (Nacalai Tesque, Kyoto, Japan) and dispersed into single cells or small clusters, which were plated onto 0.1% gelatin-coated dishes. hiPSC-CMs were maintained in a conditioned medium. The experimental study using hiPSC-CMs was approved by the Kyoto University ethics review board (G259) and conformed to the principles of the Declaration of Helsinki.
Electrophysiological recordings of hiPSC-CM APs. For single-cell patch-clamp recordings, gelatincoated glass coverslips were placed into each well of a 6-well plate. Two milliliters of DMEM/F12 containing 2% fetal bovine serum and 80,000-120,000 CMs were added to each well. Spontaneous APs were recorded from beating single CMs using the perforated patch-clamp technique with amphotericin B (

Results
Mapping the magnitude of MSE over the nine global parameter space. Parameter identifiability is a central issue in the parameter optimization of biological models 14,20 . To confirm the identifiability of a unique set of sf x s using the parameter optimization method, mapping of the MSE distribution over an enlarged parameter space defined by the sf x of the nine ionic currents of the baseline model is required. The randomization of sf x ranged from 1/10 to approximately 10 times the default values. The calculation was performed for approximately 5,000,000 sets, as shown in Fig. 1; magnitudes of log(MSE) were plotted against each sf x on the abscissa. The data points of MSE at a given sf x include all variable combinations of the other eight sf x s. The algorithm of the PS method searches for a parameter set, which gives the minimum MSE at a given stp through the process of optimization. Drawing a clear envelope curve by connecting the minimum MSEs at each sf x was difficult because of the insufficient number of data points in these graphs (Fig. 1). Nonetheless, an approximate envelope of the minimum MSEs may indicate a single global minimum of MSE located at the control sf x equals one, as typically exemplified by sf Kr -and sf bNSC -MSE relations. On both sides of the minimum, steep slopes of MSE/sf x were evident in all graphs. Outside this limited sf x -MSE area, the global envelope showed a gentle and monotonic upward slope toward the limit on the right side. No local minimum was observed in all of the sf x -MSE diagrams, except for the central sharp depression. Essentially the same finding was obtained in another hiPSC-CM model (Cell 38), which showed less negative MDP (see Fig. S4 in Supplemental Materials). It was concluded that the theoretical model of cardiac membrane excitation (hVC model) has only a single central sharp depression corresponding to the control model parameter. Figure 2 illustrates the records of spontaneous APs (red traces) obtained from 12 experiments in the maximum diastolic potential (MDP) sequence (see Supplemental Materials for details). All experimental records were superimposed with simulated AP traces (black traces) obtained using conventional manual fitting. In most cases, an MSE of 1-6 mV 2 remained (Eq. 7) at the end of the manual fitting. This extra component of MSE might be largely attributed to slow fluctuations of V m of unknown origin in experimental recordings, because the non-specific random fluctuations were quite different from the exponential gating kinetics of ion channels calculated in mathematical models. This extra noise seriously interfered with the assessment of the accuracy of the parameter optimization of G x in the present study. Thus, APs produced by the manual adjustment (cell-specific model) was used as the target AP that was completely free from the extra noise when examining the feasibility of the parameter optimization algorithm.

Necessity for parameter optimization as indicated by hiPSC-CM APs.
Comparison of the AP configurations between these hiPSC-CMs clearly indicated that the classification of these APs into atrial, ventricular, and nodal types was impractical, as has been described 7 . On the other hand, if provided with the individual models fit by objective parameter optimizing tools using the baseline model (black trace), the results should be fairly straightforward for estimating the functional expression level of ion channels and to clarify the role of each current system or the ionic mechanisms in generating the spontaneous AP configuration in a quantitative manner. Thus, the objective parameter optimization of the mathematical model is a vital requirement in cardiac electrophysiology. Table 1 lists the AP metrics, including cycle length (CL), peak potential of the plateau (OS), MDP, and AP duration measured at -20 mV in addition to the MSE between individual experimental records and the model output fitted by manual fitting. CL, MDP, and AP varied markedly among different AP recordings of cells (Fig. 2). Cells were arranged according to the MDP sequence.
Feasibility of PS algorithm for parameter optimization of membrane excitation models. Automatic parameter optimization has been applied to the model of cardiac membrane excitation in a limited number of studies (for reviews, see 23 www.nature.com/scientificreports/ models of cardiac membrane excitation composed of both ionic channel and ion transporter models, except for one pioneering study 12 , which applied a more general gradient-based optimization method to the simple ventricular cell model (Beeler and Reuter [BR] model) 40 . Figure 3 shows a typical successful run of the new PS method in hiPSC-CMs. An approximate MDP of −85 mV was evident. The PS parameter optimization was started after randomizing the sf x s of the six major currents (I Kr , I CaL , I Na , I ha , I K1 and I bNSC, ) in the manual fit model within a range of ± 15% around the default values (normalized magnitude of 1). Figure 3A 1-3-compare the simulated V m,orp (black) with the target V m,t (red) at the repeat numbers of N = 1, 50, and 1167, respectively (Eq. 8). The overshoot potential (OS), APD, and CL of spontaneous AP were markedly different during the first cycle of AP reconstruction (Fig. 3A-1 The time course of decreasing log(MSE) evoked by the multi-run PS optimization was plotted for each sf x every time the set of base points was reset (Fig. 3B-1). Figure 3B-2 shows all of the log(MSE) obtained at every    Figure 4 shows the results of the orp tests, in which the optimization shown in Fig. 3 was repeated several hundred times. All results were plotted in a common coordinate of log(MSE) and individual sf x s. The population of sf x correctly converged at a single peak point very close to 1 with increasing negativity of log(MSE) for sf Kr , sf CaL , and sf Na . In contrast, sf ha, sf K1, and sf bNSC showed obvious variance. Nevertheless, they also showed a clear trend toward convergence to 1 on average. We could find clear convergence of less number (4) of sf x s used for PO method in cells, which showed relatively low MDPs as shown in the supplemental Materials (Fig. S5). Table 2    www.nature.com/scientificreports/ magnitude from the initial level just after randomization by the orp test, as in the successful example shown in Fig. 3B. The mean of individual sf x s was very close to 1 with a minimum standard error (SE) of the mean, which was < 1% of the mean, even for I K1 , I bNSC , and I ha , which showed weak convergence against log(MSE) . These results validate the accuracy of the parameter optimization using the multi-run PS method in all 12 cell-specific models, which include many varieties of spontaneous AP recorded in hiPSC-CMs.
Complementary relationship among I K1 , I ha, and I bNSC . Figure 5A illustrates the distribution of the sf x s amplitude in the top 20 data points. The final sf x s in individual runs were connected with lines for each run of PS in Cell 86 (Fig. 2). The SEM values were quite small (< 1%) in sf Kr and sf CaL . In contrast, sf ha , sf K1 , and sf bNSC showed larger deviations. This finding is interesting because the former currents are mainly involved in   Figure 5B,C show the distribution of the sf x points in the space of the three sf x dimensions. In Fig. 5B, the 20 data points seemed to be dispersed randomly in the parameter space. However, when the space was rotated to a specific angle, a linear distribution was observed (Fig. 5C), indicating that the points were distributed approximately on a plane surface in the three-dimensional space. Using multiple regression analysis, we obtained an equation that fit the 20 data points as follows (R 2 = 0.872): By replotting the data points in the two-dimensional space with the abscissa for the sum of two inward-going currents (0.76 sf ha + 0.19 sf bNSC ) and the ordinate for the outward current 0.62 sf K1 , we obtained the regression line shown in Fig. 5D. Close correlations among the three sf x s were indicated with a high R 2 of 0.941. This finding confirms that the three currents have complementary relationships with each other to provide virtually identical configurations of spontaneous AP. In other words, log(MSE) remains nearly constant as long as the composition of the currents satisfies the relationship given by Eq. 16.
The complementary relationship was further examined by performing an orp test after fixing one of the two factors, sf K1 or (sf ha + sf bNSC ), as illustrated in Fig. 5B. Figure 6A shows the log(MSE) vs. sf K1 relation when (sf ha + sf bNSC ) was fixed at the values obtained by the orp test. Indeed, the typical convergence of the sf K1 was obtained. Alternatively, if the sf K1 was fixed, the convergence was obviously improved for both sf ha and sf bNSC (Fig. 6B-1,2), but it was less sharp if compared to sf Kr , sf CaL and sf Na (not shown, but refer to corresponding results in Fig. 4A). This finding was further explained by plotting the relationship between the two inward currents, I ha and I bNSC, as illustrated in Fig. 6C. The regression line for the data points was fitted by Eq. 17 with R 2 = 0.86, supporting the complemental relationship between the two inward currents, I ha and I bNSC .
The moderately high R 2 indicates that the SDD is determined not only by the major I ha and I bNSC but also by other currents, such as I K1 , I Kr , the delayed component of I Na (I NaL ) and I CaL , which were recorded during the SDD as demonstrated in Fig. 3.
Essentially the same results of complementary relationship among sf ha , sf bNSC and sf K1 were obtained in Cell 91, which also showed the long-lasting SDD with the very negative MDP as in Cell 86, as shown in Fig. 2 and Table 2. The regression relation for the data points was fitted by Eqs. (18) and (19) with R 2 = 0.656 and 0.472, respectively.
Principal components in the hiPSC-CM model. The PS frequently got stuck during the progress of parameter optimization and failed to reach the global minimum in the present study (Figs. 4,6). The major cause of this interruption may most probably be attributed to the fact that sf x s were used directly as the search vector of the PS. In principle, the algorithm of PS parameter optimization gives the best performance when the parameters search is conducted in orthogonal dimensions where each dimension does not affect the adjustment of other sf x 41 . To get deeper insights, we applied the principal component (PC) analysis to the set of 6 sf x s selected in (16) 0.762 · sf ha − 0.619 · sf K1 + 0.191 · sf bNSC = 0.333554 (17) 0.9736 · sf bNSC + 0.2281 · sf ha = 1.2024 (18) 0.572 · sf ha − 0.132 · sf K1 + 0.810 · sf bNSC = 1.25891 (19) 0.9279 · sf ha + 0.3706 · sf bNSC = 1.30025 www.nature.com/scientificreports/ the baseline model. We performed PC analysis on the data points recorded in the vicinity of the minima (using the top 20 data). As illustrated in Fig. 7, each of the 6 PCs was not composed of a single sf x but mostly included multiple sf x sub-components. This finding indicates the inter-parameter interactions during the process of parameter optimization. For example, the changes in sf K1 or sf bNSC simultaneously affect PCNo.1, 3, 6 or 1, 2, 3 PCs, respectively. Both sf CaL and sf Kr affect PCNo.4, 5. It might be concluded that the frequent interruptions of PS parameter optimization are most probably caused by the sporadic appearance of the local minima of MSE through interactions among sf x s.

Discussion
New findings in the present study are listed below.
(1) Mapping of the MSE distribution over the enlarged parameter space was conducted by randomizing the G x s of the baseline model. It was confirmed that the baseline model had only a single sharp depression in MSE at the default G x s (Fig. 1). (2) The preliminary cell-specific models were firstly prepared by the conventional manual tuning of G x s to superimpose the model output on each of the 12 experimental AP recordings (Fig. 2). The parameter search space was restricted to a relatively small space to facilitate parameter optimization. (3) The sf x s of the 4-6 G x parameters were initially assigned random values from a uniform distribution ranging between ± 10% of the default values. The MSE was calculated between the randomized model output and the intact model AP as the target of optimization (Fig. 3). (4) Plotting the parameter sf x in common sfx-MSE coordinates during each run of several hundred runs of optimization (Fig. 4), revealed that the sf x distributions of I Kr , I CaL , and I Na converged sharply to a single point with decreasing MSE, which exactly equaled the default values. In contrast, estimates of sf K1 , sf ha , and sf bNSC deviated slightly within a limited range around the default values in cells showing long-lasting SDD (Fig. 4). (5) For statistical evaluation, the mean ± SE of sf x in the top 20 MSE estimates was calculated for individual cells ( Table 2). The results of the parameter optimization in the 12 cells indicated that the means of sf x s were very close to 1.00, with an SE < 0.01 for all G x s. (6) A complementary relationship was found between sf K1 , sf ha , and sf bNSC in determining the gentle slope of the long-lasting SDD in two representative cells (Fig. 5). Supporting this view, sf K1 clearly focused on the unit provided that sf ha and sf bNSC were fixed and vice versa (Fig. 6). (7) The six search vectors of sf x in the presented model could be replaced by the same number of theoretical PCs, and each PC was mostly composed of multiple sf x s (Fig. 7). This finding supports the view 12 that the complex interactions among I x s might interrupt the progress of the parameter optimization when sf x s are used as the search vector instead of using theoretical orthogonal ones.
The use of an initial randomized set of parameters was crucial in examining whether an optimization method could determine unique estimates independent of the initial set of parameters, as used in the GA-based method for determining the G x s of the mathematical cardiac cell model 23 . The aforementioned seven findings confirm the feasibility of the PS method. Most likely, the PS method is applicable to variable mathematical models of other cell functions. Reference 26 provides a more systematic review of parameter optimization in cardiac model development.
It has been suggested that different combinations of parameters generate similar outputs 12,[23][24][25] . In the present study, this suggestion was explained at least in part by the complementary relationship, for example, between I K1 , I ha , and I bNSC in determining dV m /dt of SDD, which is a function of the total current (Eq. 2, Figs. 5,6). The gradient-based optimization method relies on the precise variation in the time course of dV m /dt induced by  Fig. 4, which showed good optimization results (Cell 86). Each magnitude of 6 PCs was normalized to give a unit magnitude. Note each PC is composed of multiple components of ionic current, which are indicated in the Index with corresponding colors. www.nature.com/scientificreports/ time-dependent changes in individual sf x s (Eq. 2). Therefore, MSE was calculated over the entire time course of spontaneous APs. Notably, we did not use AP metrics, which indirectly reflect the kinetic properties of the individual currents. Even with this measure of calculating the MSE, the time-dependent changes in pO (Eq. 3) might be relatively small between two major currents, I K1 and I ha , in comparison to I bNSC , which has no V m -dependent gate during the SDD, as shown in the current profile in Fig. 3A-3. We assume that the gradient-based optimization method can determine different contributions of individual currents if the optimization is conducted only within a selected time window of SDD. If the MSE is calculated over multiple phases of spontaneous AP, the influence of a particular phase on the MSE should be diluted. In our preliminary parameter optimization, this problem was partly solved using a weighted sum for different phases of the spontaneous AP in summing the MSE. The small amplitude of a given current might be an additional factor in the weak convergence of sf x observed in the diagram of sf x -MSE in the orp test of optimization. If the current amplitude is much smaller than the sum of all currents in determining dV m /dt (Eq. 2), the resolution of the PS method would decrease. Sarkar et al. 24 demonstrated that the model output, for example, the AP plateau phase, was almost superimposable when different ratios of G Kr and G pK were used to reconstruct the model output (their Fig. 1). The authors reported that the AP metrics used for comparisons, such as APD, OS, and APA, were quite similar. However, these results were obtained by applying different combinations of sf x to the same Ten Tusscher-Noble-Noble-Panfilov (TNNP) model 33 . This means that the relative amplitudes of I Kr and I pK in the TNNP model were much smaller than those of the major I CaL during the AP plateau, even though I Kr and I pK have completely different gating kinetics. Thus, the results of the parameter optimization should be model-dependent. The same arguments can also be applied to the use of the FR guinea pig model 42 in the study by Groenendaal et al. 23 .
A gradient-based parameter optimization method was applied to the cardiac model of membrane excitation in a study 12 that analyzed the classic BR model 40 . The whole-cell current in the BR model is composed of a minimum number of ionic currents, background I K1 , and three time-dependent currents (I Na , I s , and I x1 ) which were dissected from the voltage clamp experiments by applying the sucrose gap method to the multicellular preparation of ventricular tissue. The gating kinetics of the latter three currents were formulated according to Hodgkin-Huxley type gating kinetics, which is quite simple compared to the recent detailed description of ionic currents. The authors reported that the parameter optimization was difficult if the AP configuration was used as the target of the parameter optimization. They used the time course of the whole-cell current as a target for parameter optimization. However, the number of parameters was quite large (in their study, 63) and included limiting conductances and gating kinetics. The authors suggested that the feasibility of the parameter optimization method would be improved with additional experimental data.
In modern mathematical cardiac cell models, most ionic currents are identified by whole-cell voltage clamp and single-channel recordings in dissociated single myocytes 43 using the patch-clamp technique 44 and by identifying the molecular basis of membrane proteins. The molecular basis of the ion channels expressed in the hiPSC-CMs from the GSE154580 GEO Accession viewer is mostly identical to those in the adult cardiac myocytes, rather than in the fetal heart. Moreover, gating kinetics have been extensively studied to characterize ionic currents within the cell model. In principle, the detailed characterization of individual currents should facilitate the identifiability of the model parameter, but should not necessarily interfere with parameter optimization. We consider that the manual fitting of the model parameters to the AP recording using a priori knowledge of biophysical mechanisms should largely facilitate the subsequent automatic parameter optimization. Also, ionic currents left at the default values work as constraints to improve the identifiability of the target parameters.
After validating the automatic parameter optimization method, the final goal of our study was to determine the principle of ionic mechanisms that are applicable to the full range of variations in spontaneous AP records in both hiPSC-CMs and mature cardiomyocytes. The multi-run PS method was applied to the experimental AP recordings using the initial parameter sets obtained by the conventional manual fit. The protocol for measuring G x s was the same as that used in the present study, except for the use of experimental AP recordings instead of the output of the ' cell-specific model' . In our preliminary analysis, the magnitude of the individual model parameters obtained by manual tuning was corrected by < 15% by objective parameter optimization.
Finally, the ionic mechanisms underlying the SDD of variable time courses will be analyzed in a quantitative manner, for example, by using lead potential analysis 45

Limitations
There are several limitations in the present study. In general, the obvious limitations of published mathematical models of cardiac membrane excitation are caused by a shortage of functional components inherent in intact cells. For example, [ATP] i that is controlled by energy metabolism is a vital factor in maintaining the physiological function of ion channels as well as the active transport of the Na + /K + pump 46 . Moreover, most models do not account for modulation of the ion channel activity through phosphorylation of the channel proteins, detailed modulation of the channel by [Ca 2+ ] i , alterations in ion channel activity by PIP 2 47,48 , and tension of the cell membrane through changes in cell volume [49][50][51][52] . The detailed Ca 2+ dynamics of [Ca 2+ ] i are still not implemented in most cardiac cell models. These dynamics include Ca 2+ release from sarcoplasmic reticulum (SR) activated through the coupling of a few L-type Ca 2+ channels with a cluster of ryanodine receptors (RyRs) at the dyadic junction 29 and Ca 2+ diffusion influenced by the Ca 2+ -binding proteins 53 . To simulate Ca 2+ binding to troponin during the development of contraction, a dynamic model of the contracting fibers is necessary [54][55][56][57] . These limitations should be thoroughly considered when investigating pathophysiological phenomena such as arrhythmogenesis. The scope of the present study was limited to the AP configurations of hiPSC-CMs, which were assumed to be www.nature.com/scientificreports/ 'healthy' with respect to the above concerns. For example, [ATP] i , [Na + ] i , and Ca tot were kept constant, and the standard contraction model was implemented, as in the hVC model. The parameter optimization presented in this study can be achieved in a practical manner by limiting the number of unknown parameters. We determined only G x s based on the assumption that ion channel kinetics are preserved, as in hiPSC-CMs and mature myocytes. Usually, four to six ionic currents are selected for optimization. The orp method could be performed simultaneously for all nine ionic currents, as described in Eq. (1). However, the computation time was radically prolonged and the resolution was not as high as that obtained using a modest number of parameters. We consider that the determination of a limited number of G x s is relevant to solving physiological problems in terms of detailed model equations for each current system.
Although I NCX and I NaK contributed sizeable fractions of the whole-cell outward and inward currents, respectively (Fig. 3A-3), we excluded the scaling factors sf NaK and sf NCX from the parameter optimization for the sake of simplicity. Instead, the possible drift of intracellular ion concentrations was fixed during the repetitive adjustment of ionic fluxes by varying sf x , as shown in Table 2. The introduction of the empirical equations (Eq. 13 and 14) was useful for adjusting [Na + ] i and Ca tot (Table 2) so that the time course and magnitude of I NCX remained almost constant during the parameter optimization. In future studies, when the influences of varying [Na + ] i and/ or Ca tot are examined under various experimental conditions, the reference levels of [Na + ] i and/or Ca tot (stdNai and stdCatot in Eq. 13 and 14) might be replaced by experimental measurements.
For the excitation-contraction coupling and calcium-induced calcium release (CICR) in hiPSC-CMs lacking T-tubules, Koivumaki et al. 58 developed the novel Paci model of hiPSC-CM with essential features of membrane electrophysiology and intracellular CICR with the spontaneous membrane excitation (mouse fetal cell model 58 ) as a platform that can be used to facilitate the translational research from hiPSC-CMs to heart diseases 59 . This composite model demonstrated spontaneous Ca 2+ release, which occurred several tens of milliseconds before the AP, to serve as a trigger. This is different from the almost simultaneous rise in spontaneous AP and the accompanying Ca 2+ transient, as demonstrated by Spencer et al. 60 . These differences might most probably be due to the variable degrees of maturation of hiPSC-CMs used in different laboratories.
The issue of coupling CICR with cardiac membrane excitation in the absence of T-tubules has long been extensively discussed in sinoatrial (SA) node pacemaker cells. Maltsev and Lakatta proposed cell models in which APs were triggered by the gradual increase in [Ca 2+ ] in the heuristic submembrane space during the Ca 2+ -transient (the 'Ca-clock' theory) 61,62 . Himeno et al. examined this issue using patch-clamp experiments in isolated SA node pacemaker cells 63 . The authors described that the spontaneous rhythm remained intact when SR Ca 2+ dynamics were acutely disrupted by addition of high doses of a Ca 2+ -chelating agent to the cytosol. This experimental finding could be reconstructed using their SA node cell model, supporting the membrane origin of spontaneous AP generation. A more detailed and extensive theoretical study was published by Stern et al. 64 (see also Hinch et al. 4 ). The authors constructed a computational cell model that included the three-dimensional diffusion and buffering of Ca 2+ in the cytosol. The Ca 2+ -releasing couplon was located at the site of close contact of the junctional SR membrane with the cell membrane, where the individual clusters of RyRs of various sizes on the SR membrane and a few LCC on the cell membrane are functionally coupled across the nanoscale gap. Interestingly, no local Ca 2+ release occurred if the clusters of RyRs were separated by > 1 μm. However, bridging large RyR clusters to form an irregular network can lead to the generation of propagating local CICR events and partial periodicity, as observed experimentally. Considering all these experimental and theoretical findings, the issue of 'Ca-clock' is still a matter of debate. Therefore, we consider that including all these details in the hiPSC version of the hVC model is clearly beyond the scope of the present study, which aims to develop a PO method to determine the parameters of the membrane excitation model in general.
The PO method was not applied to several ionic currents in this study. For example, it was difficult to determine the kinetics of the T-type Ca 2+ channel (I CaT ; Ca V 3.1) and so it was excluded from the present study. The very fast opening and inactivation rates that have been previously described 65 suggest a complete inactivation of I CaT over the voltage range of SDD, while the sizeable magnitude of the window current that has also been described 66 suggests a much larger contribution to SDD. The kinetics of I CaT remain to be clarified through experimental examination. The sustained inward current, I st , has recently been attributed most probably to Ca V 1.3 67 , which is activated at a more negative potential range than I CaL (Ca V 1.2) 68,69 . In the present study, I bNSC was used to represent the net background conductance. However, several components of background conductance have been identified at the molecular level in mature myocytes (for a review of TRPM4, see 70 ). Experimental measurements of the current magnitude of each component are required.
Gábor and Banga indicated that the multi-run method performed well in certain cases, especially when high-quality first-order information was used, and the parameter search space was restricted to a relatively small domain 16 . Another study echoed these findings 19 . In the present study, manual fitting of the parameters (Fig. 1) was required to utilize the multi-run PS method over the restricted search space. One of the major difficulties in the manual fitting of individual G x s arose during SDD, where I Kr , I K1 , I bNSC , and I ha , in addition to I NaK and I NCX , constitute the whole-cell current (Fig. 3A-3). However, close inspection of the current components in Fig. 3A-3 provides hints on how to do with the manual fit. The transient peak of I Kr dominates the current profile during the final repolarization phase from -20 to -60 mV in all 12 hiPSC-CMs 71 , since I CaL and I Ks rapidly deactivates before repolarizing to this voltage range. I NaK and I NCX are well controlled by the extrinsic regulation in Eqs. (13) and (14). Thus, manual fitting of sf Kr was first applied to determine sf Kr . The MDP more negative than -70 mV was adjusted by the sum of time-dependent (I Kr + I K1 ) and time-independent I bNSC . Then, I Kr is deactivated when depolarization becomes obvious after the MDP, and the depolarization-dependent blocking of I K1 by intracellular substances 72 play major roles in promoting the initial linear phase of SDD. Thus, the amplitudes of sf K1 and sf bNSC may be approximated during the initial half of the SDD. The latter half of SDD, including the foot of AP (i.e., the exponential time course of depolarization toward the rapidly rising phase of AP) was mainly determined by the subthreshold V m -dependent activation of I Na (after MDP more negative than -70 mV) and/or I CaL (after MDP less www.nature.com/scientificreports/ negative than -65 mV). Thus, sf Na and sf CaL were roughly determined by fitting the foot of the AP and the timing of the rapid rising phase of AP. The plateau time course of AP is determined by sf CaL and the Ca 2+ -mediated inactivation of I CaL (parameter KL 4 ). Because the kinetics of outward currents I Kur , I Kto (endo-type), and I Ks are quite different from those of I Kr , the plateau configuration was determined incrementally by adjusting these currents. We failed to observe phase 1 rapid and transient repolarization in hiPSC-CMs (Fig. 2), which is a typical sign of the absence of epicardial-type I Kto .
In hiPSC-CMs showing less negative MDP than approximately -65 mV, the contribution of I K1 , I Na , and I ha should be negligibly small because I K1 is nearly completely blocked by intracellular Mg 2+ and polyamines, I Na is inactivated, and I ha is deactivated during SDD, even if it is expressed.
Nevertheless, parameter optimization might be laborious and time-consuming for those unfamiliar with the electrophysiology of cardiac myocytes. This difficulty might be largely eased by accumulating both AP configurations and the underlying current profile obtained in parameter optimization into a database in the future. If this database becomes available, computational searches for several candidate APs for the initial parameter set will be feasible, which will be used for automatic parameter optimization.

Data availability
The AP records used in Section 4.2, and the source code of the optimization program are available in the Supplemental Material link for the following bioRxiv entry. https:// doi. org/ 10. 1101/ 2022. 05. 16