Introduction {#Sec1}
============

Several hundred types of cells develop from a single genome through accurate spatiotemporal regulation of gene expression. The vertebrate heart is a good example of this phenomenon, as it substantially changes its shape and function at the cellular, tissue, and organ levels throughout a lifetime. In early embryonic development, the heart becomes a functional organ, acting as a pump. The heart develops and gains new functions while continuously pumping blood, and heart abnormalities during the early developmental stages progress to congenital heart malformations; therefore, the developmental program of the heart, including the expression of the genes responsible for various ionic channels, is likely to be tightly regulated. Electrophysiological recordings of various ionic channels and quantification of the genes responsible for the channels have been reported primarily for 4 representative stages: early embryonic (EE), late embryonic (LE), neonatal, and adult. To provide a complete overview of developmental regulation, it is necessary to observe the developmental changes occurring in the heart across these representative stages.

In rodents, spontaneous action potentials (APs) have been reported for the EE stage in developing rodent ventricular myocytes, eventually disappearing in passive contracting cells in the LE stage \[[@CR1]\]. The electrophysiological properties of individual ion channels have been investigated in isolated ventricular myocytes at the 4 representative stages by means of patch-clamp methods \[[@CR2]--[@CR4]\]. In addition to cells that exhibit spontaneous APs, embryonic rodent ventricular tissues also contain quiescent cells with no spontaneous APs. In 12-day fetal hearts, for example, 9 of 14 ventricular cells were quiescent and exhibited a resting membrane potential (RMP) of −48.4 ± 1.8 mV \[[@CR5]\]; similarly, in 18-day postcoitum (dpc) mice, 6 of 13 isolated ventricular cells were spontaneously active, whereas the remaining 7 cells were quiescent \[[@CR6]\]. Moreover, the beating rate of the spontaneous APs ranged from 35 ± 11 beats per minute (bpm) \[[@CR5]\] to 232 bpm \[[@CR7]\] in 12.5-dpc embryonic rat ventricles, and from 178 ± 12.7 \[[@CR8]\] to 124 ± 8.7 bpm \[[@CR9]\] in 9.5-dpc embryonic mouse ventricles; the reported beating rates roughly correspond to a basic cycle length (BCL) of approximately 259--2,500 ms. In addition to regular spontaneous APs, irregular spontaneous APs have also been reported in the embryonic ventricular cells of both mice \[[@CR6]\] and rats \[[@CR5]\]. Previously, we modeled the developmental changes in the APs of cardiac ventricular myocytes \[[@CR10]\] by using the Kyoto model \[[@CR11]\] and the Luo--Rudy dynamic (LRd) model \[[@CR12]\]. The measured APs at developmental stages were reproduced using common sets of these models by varying the relative densities of the ionic currents, pumps, exchangers, and sarcoplasmic reticulum (SR) Ca^2+^ kinetics. In addition, Jonsson et al. \[[@CR13]\] combined molecular biology and computer simulations to demonstrate that human embryonic stem cell-derived cardiomyocytes (hESC-CMs) have an immature electrophysiological phenotype, based on insufficient function of inward rectifier K^+^ current (*I*~K1~) channels and a shift in the activation of sodium channels. Thus, computer simulation is a powerful approach for confirming experimental data and providing insights into the possible functional mechanisms involved in cardiac development. Although our simulations well reproduced the disappearance of the spontaneous APs between the EE and LE stages as well as the changes in AP durations during the transition from the LE to the neonatal and adult stages, we have not clarified the contribution of each ionic current system to the reported characteristics of rodent ventricular cells.

In the present study, we examined the functional changes in developing embryonic ventricular cells to identify the pivotal components in the model in order to describe the reported characteristics of embryonic rodent ventricular cells. We switched the relative densities of ionic components that differ between the EE and LE stages, and tested 512 combinations with the Kyoto model, 128 combinations with the Ten Tusscher--Panfilov (TP) human ventricular cell model \[[@CR14]\], and 32 combinations with the LRd model. The 160 regular spontaneous APs predicted in the Kyoto model had a wide range of BCL values, all of which were within the reported range of the BCL in 9.5-dpc mice \[[@CR6], [@CR8], [@CR9]\] and in 11.5-dpc rats \[[@CR5], [@CR7]\]. In all 3 models, the combinations in which *I*~K1~ was increased before the disappearance of *I*~f~ were predicted to result in high \[Ca^2+^\]~i~, and most of the combinations had quiescent membrane potentials slightly positive to −80 mV, as reported in 12.5-dpc fetal rat ventricular cells \[[@CR5]\].

Methods {#Sec2}
=======

Previously, we simulated the APs of rodent ventricular cells at the EE, LE, and neonatal stages by using the Kyoto model---an electrophysiological model of guinea pig ventricular cells \[[@CR10], [@CR11]\]. Briefly, quantitative changes in various ionic components were represented as the densities of the components in the developmental stages relative to those in the adult stage. These relative densities were then multiplied by the corresponding conductance (nS/pF) or conversion factors (pA/pF·mM) to demonstrate that developmental changes in the APs can be reproduced using common sets of mathematical equations. We adopted the same procedure in the present study, and reconstructed EE and LE ventricular cell models by using the updated Kyoto model \[[@CR15]\].

Ionic currents, exchangers, and SR Ca^2+^ kinetics {#Sec3}
--------------------------------------------------

Quantitative changes in ionic currents (Fig. [1](#Fig1){ref-type="fig"}) were either computed from current--voltage (*I*--*V*) curves or estimated on the basis of qualitative observations as relative densities. In our previous study \[[@CR10]\], we showed that the EE stage corresponds to approximately 9.5 dpc in mice and 11.5 dpc in rats, and the LE stage corresponds to 1--5 days before birth; we used the term "rodent" to represent all species from which we obtained experimental data; however, this does not imply that the models represent rodent ventricular cells in general. Table [1](#Tab1){ref-type="table"} lists the membrane currents that quantitatively change during embryonic development. The relative densities were as reported previously \[[@CR10]\], except for the rapid component of the delayed rectifier K^+^ current (*I*~Kr~) and the background nonselective cation current (*I*~bNSC~). The relative density of *I*~Kr~ in the EE ventricular cell was originally set to 10.0 times that of the adult stage in our previous study \[[@CR10]\], on the basis of the qualitative observation using the selective *I*~Kr~ blocker in EE rats \[[@CR16], [@CR17]\], and the models were constructed on the basis of the previous version of the Kyoto model \[[@CR11]\]. In this study, we set the relative density of *I*~Kr~ in the EE stage to 2.0 times that of the adult stage because the amplitude of *I*~Kr~ conductance was increased by fourfold to reconstruct the AP duration prolongation through the pharmacological inhibition of *I*~Kr~ in the current version of the Kyoto model \[[@CR15]\]. The relative density of *I*~Kr~ in the LE stage was also set to 2.0 times that of the adult stage on the basis of the qualitative observations using the selective *I*~Kr~ blocker in LE rats \[[@CR18]\] and the *I*--*V* curve of the sum of *I*~Kr~ and *I*~Ks~ in LE guinea pigs \[[@CR2]\]. *I*~bNSC~ was assumed to have a constant current density during development.Fig. 1Developmental changes in ionic components of rodent ventricular cells. Relative densities of the Na^+^ current (*I* ~Na~), L-type Ca^2+^ current (*I* ~CaL~), inward rectifier K^+^ current (*I* ~K1~), transient outward current (*I* ~to~), and ATP-sensitive K^+^ current (*I* ~KATP~) increase during embryonic development (see Table [1](#Tab1){ref-type="table"} for details). The funny current (*I* ~f~) and sustained inward current (*I* ~st~) were assumed to completely disappear in late embryonic (LE) ventricular cells (Table [2](#Tab2){ref-type="table"}). The relative densities of the Na^+^/Ca^2+^ exchange current (*I* ~NaCa~) and sarcoplasmic reticulum (SR) components are known to change reciprocally \[[@CR20]\]; the relative density of *I* ~NaCa~ decreases as the SR-related proteins, ryanodine receptor (RyR) proteins and SERCA, increase (Table [3](#Tab3){ref-type="table"}). The *gray bars* indicate the gestational durations of rat \[[@CR34]\], rabbit \[[@CR21]\], and guinea pig \[[@CR35]\]. We defined that the EE stage corresponds to approximately 9.5 dpc in mice and 11.5 dpc in rats, and the LE stage corresponds to 1--5 days before birthTable 1Relative densities for ionic currents obtained from the literatureEEReferenceLEReference*I* ~Na~0.07\[[@CR3]\]1.00\[[@CR3]\]*I* ~CaL~0.46\[[@CR22]\]0.78\[[@CR2]\]*I* ~K1~0.11\[[@CR36]\]1.00\[[@CR2]\]*I* ~Kr~2.0\[[@CR16], [@CR17]\]2.0\[[@CR2], [@CR18]\]*I* ~KATP~0.32\[[@CR3], [@CR37]\]0.88\[[@CR37], [@CR38]\]*I* ~to~0.01\[[@CR3]\]0.27\[[@CR38]\]The densities relative to that in the adult stage of the Na^+^ current (*I* ~Na~), L-type Ca^2+^ current (*I* ~CaL~), inward rectifier K^+^ current (*I* ~K1~), ATP-sensitive K^+^ current (*I* ~KATP~), and transient outward current (*I* ~to~) for the early embryonic (EE) and late embryonic (LE) stages were estimated from the current--voltage (*I*--*V*) curves of cells in vitro. The *I*--*V* curve for *I* ~Na~ was obtained from 11- to 13-day postcoitum (dpc), and 17- to 20-dpc mice \[[@CR3]\]. For *I* ~CaL~, the early embryonic *I*--*V* curve was obtained from 9.5-dpc mice \[[@CR22]\]; the late embryonic *I*--*V* curve was obtained from 18-dpc mice \[[@CR22]\] and fetal guinea pigs 1--7 days before birth \[[@CR2]\]. For *I* ~K1~, *I*--*V* curves were obtained for the 12-dpc rat \[[@CR36]\] and the fetal guinea pig 1--7 days before birth \[[@CR2]\]. The relative densities of *I* ~KATP~ were obtained on the basis of the data for the 12-dpc and 18-dpc rats \[[@CR37]\]. For *I* ~to~, *I*--*V* curves were obtained for 11-dpc mice \[[@CR3]\] and 1-day-old rats \[[@CR38]\]

As listed in Table [2](#Tab2){ref-type="table"}, the conversion factors for *I*~f~ (*P*~f,K~ and *P*~f,Na~) were set to those of *I*~f~ in sinoatrial node (SAN) cells \[[@CR19]\] for the EE stage and to 0 pA/pF mM for the LE stage, based on the *I*--*V* curves for *I*~f~ in embryonic mice \[[@CR6]\]. The conversion factors for *I*~st~ (*P*~st,K~ and *P*~st,Na~) were arbitrarily set to half of those in SAN cells for the EE stage and to 0 pA/pF mM for the LE stage because there was no evidence for the presence of *I*~st~ in embryonic rodents.Table 2Estimated conversion factors for the funny current (*I* ~f~) and sustained inward current (*I* ~st~)SAN\
(pA/pF mM)ReferenceEE\
(pA/pF mM)ReferenceLE\
(pA/pF mM)Reference*P* ~f,Na~0.01379\[[@CR19]\]0.01379\[[@CR6]\]0\[[@CR6]\]*P* ~f,K~0.05635\[[@CR19]\]0.05635\[[@CR6]\]0\[[@CR6]\]*P* ~st,Na~0.007375\[[@CR19]\]0.003588Arbitrarily set0Arbitrarily set*P* ~st,K~0.0043125\[[@CR19]\]0.002156Arbitrarily set0Arbitrarily setThe conversion factors for *I* ~f~ and *I* ~st~ and the permeability coefficients for Na^+^ and K^+^ (*P* ~f,Na~, *P* ~f,K~, *P* ~st,Na~, *P* ~st~ ~K~) for the EE and LE stages were arbitrarily set on the basis of the *P* ~f,Na~ and *P* ~f,K~ in the sinoatrial node (SAN) model \[[@CR19]\] and the expression of the mRNA encoding the *I* ~f~ channels in the 9.5-dpc and 18-dpc mice \[[@CR6]\]. The *P* ~st,Na~ and *P* ~st,K~ in the EE model were assumed to be approximately half those in the SAN model \[[@CR19]\]

SR development in ventricular cells of developing rodents is highly correlated with gestational age (Fig. [1](#Fig1){ref-type="fig"}): guinea pigs, which have a fairly long gestation period (67 days), have an almost fully developed SR at birth, whereas rabbits and rats have an immature SR that increases in size during postnatal development by 3- and 5-fold, respectively, and reaches adult size in a few weeks \[[@CR20]\]. The relative amount of SR Ca^2+^ pump proteins in rabbit \[[@CR21]\] and mouse \[[@CR22]\] and that of the ryanodine receptor (RyR) channel in mouse \[[@CR22]\] were adopted to be representative of SR Ca^2+^ kinetics for the rodent LE stage (Table [3](#Tab3){ref-type="table"}). Similarly, we computed the average relative amount of the Na^+^/Ca^2+^ exchanger in rabbit \[[@CR20]\] and mouse \[[@CR22]\]. In addition, the Ca^2+^-induced activation rate for the RyR channel was multiplied by the average relative density values for SR Ca^2+^ kinetics to represent changes in the Ca^2+^-induced Ca^2+^ release (CICR) factor during embryonic development \[[@CR4]\].Table 3Relative ratios of ion fluxes in exchanger, pump, and sarcoplasmic reticulum (SR) Ca^2+^ kineticsEEReferenceLEReferenceNa^+^/Ca^2+^ exchange4.95\[[@CR22]\]2.0\[[@CR20], [@CR22]\]SR Ca^2+^ pump0.03\[[@CR22]\]0.21\[[@CR21], [@CR22]\]RyR channel0.05\[[@CR22]\]0.40\[[@CR22]\]SR Ca^2+^ transfer0.04Arbitrarily set0.30Arbitrarily setSR Ca^2+^ leak0.04Arbitrarily set0.30Arbitrarily setCICR factor−3Arbitrarily set−60Arbitrarily setDevelopmental changes in Na^+^/Ca^2+^ exchange (*I* ~NaCa~) have been reported as western blots of the NCX1 protein in rabbit \[[@CR20]\] and mouse \[[@CR22]\]. Based on findings implying that postnatal quantitative changes in the density of *I* ~NaCa~ are in good agreement with changes in protein production levels \[[@CR35]\], we assumed that the relative amounts of the proteins directly reflected the relative ratios of the ion fluxes of the *I* ~NaCa~, SR Ca^2+^ pump, and RyR channel. Hence, the western blots of the SR Ca^2+^ pump proteins in rabbit \[[@CR21]\] and mouse \[[@CR22]\] and that of the RyR channel in mouse \[[@CR22]\] were adopted for quantifying the relative amount of ion fluxes through the SR Ca^2+^ pump and RyR channel, respectively. Similarly, we computed the relative amount of ion flux through Na^+^/Ca^2+^ exchange on the basis of the relative amount of the Na^+^/Ca^2+^ exchanger protein in rabbit \[[@CR20]\] and mouse \[[@CR22]\]. The relative fluxes of Ca^2+^ transfer from the SR uptake site to the release site and Ca^2+^ leak from the SR were both set to 0.04 for the EE stage and 0.30 for the LE stage, based on the average of the relative amount of the RyR channel and that of the SR Ca^2+^ pump at the corresponding stages. Levels of the Ca^2+^-induced Ca^2+^ release (CICR) factor in the EE and LE stages were determined on the basis of the average relative density values for SR Ca^2+^ kinetics \[[@CR4]\]

Volumes of cell compartments {#Sec4}
----------------------------

In the Kyoto model, 80 % of the total cell volume (*V*~t~, 16,000 μm^3^) is considered accessible for ion diffusion (*V*~i~, 12,800 μm^3^) \[[@CR15]\]. The volumes of the cell compartments of the EE ventricular cell were computed as described previously \[[@CR10]\]; the volumes of the SR uptake site (*V*~up~) and the SR release site (*V*~rel~) were set to 3.394 and 1.3576 μm^3^, respectively, based on the *V*~i~ in the EE ventricular cell model (*V*~i,EE~, 1697 μm^3^). In addition, the updated Kyoto model included the mitochondrion as a cell compartment in addition to the SR, and the mitochondrial volume (*V*~mit~) was set to 23 % of *V*~i~ \[[@CR15]\]. Because the mitochondria are marginally developed in rodent EE ventricular cells, we arbitrary set the *V*~mit~ for the EE ventricular cells to 2.3 % of *V*~i,EE~ (39.031 μm^3^), which is 1/10 the ratio in the adult stage.

Switching the relative densities of the ionic components {#Sec5}
--------------------------------------------------------

We selected the following nine components to be switched between the EE and LE stages: *I*~f~, *I*~st~, *I*~K1~, Na^+^ current (*I*~Na~), L-type Ca^2+^ current (*I*~CaL~), Na^+^/Ca^2+^ exchange current (*I*~NaCa~), transient outward current (*I*~to~), ATP-sensitive K^+^ current (*I*~KATP~), and a set of 4 electrical components of the SR. The electrical components of the SR---which included the permeability of Ca^2+^ release from SR to the dyadic space through the RyR channel (*I*~RyR~), Ca^2+^ leak from the SR (*I*~SR,leak~), the SR Ca^2+^ pump, and Ca^2+^ transfer from the SR uptake site to the release site (*I*~SR,transfer~)---were treated as a set of components in the SR because all 4 components are located in the SR and develop along with the development of the SR. The other ionic components in the model were assumed to have constant current densities during embryonic development.

We applied the exact same procedure to simulations with the TP human ventricular cell model \[[@CR14]\], in which the following 7 components were switched between the EE and LE stages: *I*~Na~, *I*~f~, *I*~K1~, *I*~CaL~*I*~st~, *I*~NaCa~, and SR-related components. Although the original TP model does not contain *I*~f~, we implemented a mathematical model for *I*~f~ \[[@CR23]\]. Similarly, we implemented the same mathematical model for *I*~f~ \[[@CR23]\] in the LRd model \[[@CR12]\] and switched the *I*~Na~, *I*~f~, *I*~K1~, *I*~CaL~, and SR-related components between the EE and LE stages for further confirmation of our simulation with the Kyoto and the TP models.

Model assumptions {#Sec6}
-----------------

We assumed that the 9 components in the Kyoto model switched relative densities directly from EE to LE values without intermediate levels, independently from the other components; therefore, we first tested 512 (2^9^) combinations of the model. We then classified the simulation results for the 512 combinations according to their electrical activities and also compared the simulated results in terms of contractile force. The Kyoto model adopted a 4-state contraction model \[[@CR24]\] to simulate cardiac cell contraction; the authors of the Kyoto model also assumed that all transition steps from cross-bridge-formed states (\[T\*\] and \[TCa\*\]) to cross-bridge-released states (\[T\] and \[TCa\]) are ATP-dependent, because ATP binding to a myosin head disrupts the cross-bridge formation between myosin and actin \[[@CR24], [@CR25]\]. Although we adopted the value of half sarcomere length (hSL, μm)---computed using the model proposed by Negroni and Lascano \[[@CR23]\]---as a quantitative parameter for cell contraction, this value is not completely quantitatively accurate because we did not consider the developmental changes in contractile proteins.

The Kyoto model contains a β1-adrenergic signaling cascade in which binding of isoproterenol (Iso) to the β1-adrenergic receptor activates protein kinase A (PKA) through activation of adenylate cyclase, and the activated PKA modulates *I*~CaL~, *I*~Ks~, phospholamban, SR Ca^2+^-ATPase (SERCA), and plasma membrane Ca^2+^-ATPase \[[@CR15]\]. The application of iso-enhanced *I*~CaL~ density had a minimal effect on the EE ventricular myocytes but had strong effects on LE and adult ventricular myocytes, as observed in an experimental study \[[@CR26]\]. Differences in β1-adrenergic modulation between EE and LE ventricular myocytes were not considered in this study because we focused on membrane excitation and contraction to present an overview of the functional landscape of developmental changes in embryonic ventricular cells.

Computer simulation procedures {#Sec7}
------------------------------

We switched the relative conductances of the 9 components between EE and LE values, and simulated the 512 combinations of the Kyoto model for 600 s with no external stimuli. For those combinations that showed no spontaneous activity for 600 s, we applied an external stimulation of −38 pA/pF, which is −8,000 pA divided by cell capacitance of the original Kyoto model (211.2 pF), at 2.5 Hz to determine whether the cells functioned as passive contracting cells. The exact same simulation procedures were applied to the TP and LRd models; we first simulated 128 combinations of the TP model and 32 combinations of the LRd model for 600 s, and provided additional 600-s simulations with external stimulation for those combinations without spontaneous activities. The external stimulation of −80 pA/pF was applied at a frequency of 2.5 Hz for the LRd model and that of −52 pA/pF at 1.0 Hz was applied for the TP model. The amplitudes of the external stimulation differed among the 3 models because we adopted the amplitudes that were used to evoke APs in the original models \[[@CR12], [@CR14], [@CR15]\].

Results {#Sec8}
=======

Classification of the 512 combinations simulated with the Kyoto model {#Sec9}
---------------------------------------------------------------------

Of the 512 combinations simulated using the Kyoto model, 248 combinations were predicted to result in quiescent cells with no spontaneous activity; the external stimulus was applied at a frequency of 2.5 Hz for 600 s to pace the 248 combinations; however, 32 of them failed to fire APs. The evoked activities are illustrated as blue hysteresis loops in Suppl. Figs. 1 and 2; the loops begin at the upper left representing their RMP, and the membrane is then depolarized, overshooting the potential, and gradually repolarized to the RMP. The sarcomere length shortens during the repolarization phase; thus, the height of the hysteresis loop represents the force of contraction. In 64 combinations in which the relative densities of *I*~Na~, *I*~f~, and *I*~K1~ were set to LE values, an increase in the relative activities of *I*~CaL~ and SR-related components resulted in larger amplitudes of hSL and a decrease in *I*~NaCa~ resulted in smaller amplitudes of hSL. The exact same results were observed in the TP model; the amplitudes of Ca^2+^ transients were increased as the relative densities of *I*~CaL~ and SR-related components were increased to the LE values, and the amplitudes were decreased as the relative *I*~NaCa~ density was decreased to the LE values (Suppl. Fig. 2).

We observed 160 regular spontaneous APs, 96 spontaneous oscillations with long AP duration, and 8 burst-like APs. The BCL of the 160 regular spontaneous APs ranged from 306 to 884 ms. The spontaneous activities are illustrated as red hysteresis loops in Suppl. Figs. 1 and 2. Because the values of the *I*~KATP~, *I*~to~, *I*~NaCa~, and SR-related components did not significantly influence the simulation results in terms of spontaneous electrical activities, we focused on the remaining 5 pivotal currents---*I*~Na~, *I*~f~, *I*~K1~, *I*~st~, and *I*~CaL~---and 32 combinations (Fig. [2](#Fig2){ref-type="fig"}) as representatives of the 512 combinations.Fig. 2Simulated membrane potential and half sarcomere length (hSL) of 32 representative combinations in which the density values for 5 components were switched. The relative densities of *I* ~Na~, *I* ~f~, *I* ~K1~, *I* ~st~, and *I* ~CaL~ were switched independently, yielding 32 combinations. The relative densities of the remaining 4 components were fixed to their EE values. Membrane potentials are presented on the *horizontal axis* and half sarcomere lengths are presented on the *vertical axis*. Models with membrane potential activities are presented as hysteresis loops, beginning at the *upper left* representing their maximum diastolic potential (MDP) and turning clockwise. Burst-like action potentials (APs) are highlighted with *dashed boxes*, shown in detail in Fig. [3](#Fig3){ref-type="fig"}a. Similarly, the combinations that were predicted to result in the spontaneous oscillation of the membrane potential with long AP duration are highlighted with the *solid black boxes*, shown in detail in Fig. [4](#Fig4){ref-type="fig"}c. The *red*hysteresis loops are illustrated on the basis of profiles from 600 to 605 s for the combinations with the regular spontaneous APs; the burst-like APs and the spontaneous oscillations with long AP duration were illustrated on the basis of the additional 600-s simulations. The *blue* hysteresis loops are illustrated on the basis of profiles from 599 to 600 s in the additional simulation with an external stimulus of 38 pA/pF at 2.5 Hz for 600 s

The changes in *I*~f~ and *I*~K1~ in other models---the TP and LRd models---demonstrated similar results in terms of regular spontaneous APs (Suppl. Figs. 3 and 4); the shifts in relative densities of *I*~f~ and *I*~K1~ to LE values terminated the spontaneous APs, which were observed when the relative densities of both *I*~f~ and *I*~K1~ were set to the EE values. The APs were not evoked by the external stimulus, when the relative densities of *I*~K1~ and *I*~st~ were set to the LE values and those of *I*~f~ and *I*~CaL~ were set to the EE values in the Kyoto model (Fig. [2](#Fig2){ref-type="fig"}, solid light boxes); similar results were observed in both the TP and LRd models, when the relative density of *I*~K1~ was set to the LE value and that of *I*~f~ was set to the EE value, regardless of the remaining components (Suppl. Figs. 3 and 4). The spontaneous oscillations with long AP duration were observed only in the Kyoto model (Fig. [2](#Fig2){ref-type="fig"}, solid black box). The burst-like behaviors were observed in 10 combinations simulated using the TP model; in the LRd model, however, all 16 combinations for which the relative density of *I*~K1~ was set to the EE value showed regular spontaneous APs, and the burst-like behaviors were not observed.

Burst-like membrane potentials in the Kyoto and TP models {#Sec10}
---------------------------------------------------------

Simulations for eight combinations demonstrated burst-like APs when the relative densities of *I*~st~ and *I*~CaL~ were set to the LE values and those of *I*~Na~, *I*~f~, *I*~K1~, and *I*~NaCa~ were set to the EE values in the Kyoto model (Suppl. Fig. 1). In the representative combination in Fig. [3](#Fig3){ref-type="fig"}a, the duration of repetitive bursts was approximately 30 s and the amplitude of the membrane potential was approximately 100 mV. The interval between the bursts was approximately 70 s at −50 mV. Similar burst-like APs were also observed in simulations using the TP model (Fig. [3](#Fig3){ref-type="fig"}b), when the relative densities of *I*~f~ and *I*~NaCa~ were set to the LE values and those of *I*~K1~ and *I*~st~ were set to the EE values (Suppl. Fig. 3); the duration of repetitive bursts was approximately 20 s and the amplitude of the membrane potential was approximately 90 mV. The interval between the bursts was approximately 80 s at −70 mV.Fig. 3Burst-like membrane potentials in the Kyoto and Ten Tusscher--Panfilov (TP) models. **a** Burst-like activity was observed when the relative densities of *I* ~st~ and *I* ~CaL~ were set to the LE values and those of the other components were set to the EE values in the Kyoto model. **b** The burst-like activity was observed when the relative densities of *I* ~f~ and *I* ~NaCa~ were set to the LE values and those of *I* ~K1~ and *I* ~st~ were set to the EE values in the TP model

The \[Na^+^\]~i~ was increased during repetitive bursts and decreased during the quiescent state between bursts in both the Kyoto and the TP models. In the Kyoto model, the \[Na^+^\]~i~ was 4.14 mM when the relative densities of all ionic components were set to EE values; the \[Na^+^\]~i~ was 3.2 mM when the burst was terminated and decreased to 2.7 mM during the quiescent state. In the TP model, however, the \[Na^+^\]~i~ was 4.01 mM when the burst was terminated and decreased to 3.65 mM, which are both higher than 1.56 mM---the \[Na^+^\]~i~ when all relative densities were set to EE values. In both models, however, the amplitudes of *I*~NaK~ increased during repetitive bursts and gradually decreased during the quiescent state between bursts.

Increase in *I*~K1~ before the disappearance of *I*~f~ resulted in high intracellular Ca^2+^ concentrations in all 3 models {#Sec11}
---------------------------------------------------------------------------------------------------------------------------

In the Kyoto model, increase in *I*~K1~ to the LE value before the disappearance of *I*~f~ resulted in either the spontaneous oscillation of the membrane potential with long APs or quiescent membrane potentials at approximately −50 mV. Of the 128 combinations in which *I*~K1~ was increased before the disappearance of *I*~f~, 32 combinations in which the relative density of *I*~CaL~ was set to the EE value and that of *I*~st~ was set to the LE value were predicted to result in quiescent membrane potentials at approximately −50 mV, and the APs were not evoked upon external stimulus application (Fig. [2](#Fig2){ref-type="fig"}, solid light boxes). Similarly, in the TP and LRd models, the APs were not evoked upon external stimulus application when *I*~K1~ was increased before the disappearance of *I*~f~ (Suppl. Figs. 3 and 4). The average RMPs were approximately −62.9 mV in the TP model and −69.1 mV in the LRd model.

In all 3 models, increase in *I*~K1~ to the LE value before the disappearance of *I*~f~ resulted in high \[Ca^2+^\]~i~ (Fig. [4](#Fig4){ref-type="fig"}a). The \[Ca^2+^\]~i~ was 12.1 μM in the LRd model and 6.27 μM in the TP model, and the decrease in the relative *I*~f~ densities decreased the \[Ca^2+^\]~i~ to 0.078 and 0.0065 μM, respectively. In the Kyoto model, we computed the average \[Ca^2+^\]~i~ during the additional 600-s simulations because \[Ca^2+^\]~i~ was not constant owing to the spontaneous oscillation of the membrane. The average of the simulated \[Ca^2+^\]~i~ was 2.38 μM when the relative *I*~f~ density was set to the EE value, and the average \[Ca^2+^\]~i~ increased to 3.48 μM as the relative *I*~f~ density decreased to 0.3; the \[Ca^2+^\]~i~ was then decreased to 0.0053 μM as the relative *I*~f~ density decreased from 0.3 to 0 (Fig. [4](#Fig4){ref-type="fig"}a). The contribution of the developmental changes in *I*~f~ and *I*~K1~ to \[Ca^2+^\]~i~ in the Kyoto model are further demonstrated in Fig. [4](#Fig4){ref-type="fig"}b, in which we shifted the relative densities of *I*~f~ and *I*~K1~ from EE to LE values independently by a 10 % increment. The average \[Ca^2+^\]~i~ was \<0.05 μM when the relative *I*~K1~ density was set to \<0.29 (30 % shift to LE), regardless of the relative *I*~f~ density; the \[Ca^2+^\]~i~ gradually increased from 0.35 to 2.38 μM as the relative density of *I*~K1~ was increased from 0.29 to 1.0.Fig. 4High intracellular Ca^2+^ concentrations were induced by an increase in *I* ~K1~ before the disappearance of *I* ~f~ in all 3 models. **a** The relative *I* ~K1~ density was set to 1.0 (LE value) and the relative densities for all the other ionic components were set to the EE values. The relative *I* ~f~ density was shifted by a 10 % interval from 1.0 to 0. For the Kyoto model, we calculated the average \[Ca^2+^\]~i~ during the additional 600-s simulation because \[Ca^2+^\]~i~ was not constant owing to the spontaneous oscillation of the membrane potential. **b** The relative densities of *I* ~f~ and *I* ~K1~ were shifted from EE to LE values independently by a 10 % interval, and the simulated \[Ca^2+^\]~i~ values are represented by *color* (range of *color bar*: 0.005--3.484 μM). **c** The spontaneous oscillation of membrane potential with long action potential duration was observed when the relative density of *I* ~K1~ was set to the LE value and those of the other components were set to the EE values. The membrane potential oscillated between approximately −35 and −80 mV

In the Kyoto model, the increase in \[Ca^2+^\]~i~ resulted in the spontaneous oscillations of the membrane potential with long APs (Fig. [4](#Fig4){ref-type="fig"}c); the membrane potential spontaneously oscillated approximately every 20 s between −35 and −80 mV. The \[Na^+^\]~i~ oscillated between 24.2 and 24.9 mM, which is considerably higher than the original \[Na^+^\]~i~ (4.14 mM) in the original EE model; the \[Na^+^\]~i~ was also increased from 9.16 to 68.9 mM in the LRd model and from 1.56 to 60.9 mM in the TP model (data not shown). The \[Ca^2+^\]~i~ also oscillated at high concentrations, from 0.46 to 3.9 μM. The membrane was depolarized to −35 mV when *I*~f~ began to apply an inward current, followed by activation of the Ca^2+^-activated background cation current (*I*~LCCa~), which is activated when the intracellular Ca^2+^ concentration (\[Ca^2+^\]~i~) is high. The membrane potential was maintained at −35 mV for approximately 20 s when the amount of the outward K^+^ current, i.e., the sum of *I*~Kr~ and *I*~K1~, was approximately the same as that of *I*~LCCa~. The rapid increase in *I*~K1~ followed by deactivation of *I*~LCCa~ subsequently led to the repolarization of the membrane to −80 mV.

Contribution of *I*~Na~ and *I*~f~ to the BCL of regular spontaneous APs {#Sec12}
------------------------------------------------------------------------

Of the 512 combinations simulated using the Kyoto model, 160 combinations were predicted to result in regular spontaneous APs. The BCL of the 160 regular spontaneous APs ranged from 306 to 884 ms; the maximum diastolic potential (MDP) ranged only from −85.0 to −80.9 mV, and the overshoot potential ranged from 1.6 to 54.3 mV. Of the 9 components that were shifted between EE and LE values, the developmental changes in *I*~Na~ and *I*~f~ made large contributions to the variation in the BCL and the overshoot potential in the Kyoto model; therefore, we further tested the contribution of the developmental changes in *I*~Na~ and *I*~f~ to the BCL and the overshoot potential.

As shown in Fig. [5](#Fig5){ref-type="fig"}a, the BCL of the regular spontaneous APs ranged from 314 to 1,550 ms in 121 combinations for which relative densities of *I*~Na~ and *I*~f~ were shifted independently by a 10 % increment from EE to LE values. The spontaneous APs ceased when the relative *I*~f~ density was set to 0 (LE value) and the relative *I*~Na~ density was set to \<0.44 (40 % shift to LE), which are illustrated as gray boxes in Fig. [5](#Fig5){ref-type="fig"}. The increase in relative *I*~Na~ density from 0.07 to 1.0 shortened the BCL, and the decrease in relative *I*~f~ density from 1.0 to 0 prolonged the BCL; therefore, the longest BCL (1550 ms) was observed when the relative *I*~Na~ density was set to 0.07 (EE value) and the relative *I*~f~ density was set to 0.1 (90 % shift to LE). The overshoot potential, on the other hand, ranged from 2.1 to 53.7 mV in the 121 combinations, and both the increase in relative *I*~Na~ density and the decrease in relative *I*~f~ density increased the overshoot potential (Fig. [5](#Fig5){ref-type="fig"}b).Fig. 5Contributions of *I* ~Na~ and *I* ~f~ to basic cycle length and overshoot potential of regular spontaneous action potential. The relative densities of *I* ~f~ and *I* ~K1~ were shifted from EE to LE values independently by a 10 % interval, and the simulated basic cycle lengths (**a**) and overshoot potential (**b**) were represented by *color*; range of color bar: (**a**) 300--1,200 ms and (**b**) 0--53.97 mV. The *gray boxes* represent the combinations with which spontaneous action potentials were terminated

Representative developmental changes in APs as *I*~Na~, *I*~f~, and *I*~K1~ were sequentially switched {#Sec13}
------------------------------------------------------------------------------------------------------

Figure [6](#Fig6){ref-type="fig"} summarizes the changes in the APs and accompanying ionic currents as we sequentially switched *I*~Na~, *I*~f~, and *I*~K1~ from the EE model; spontaneous APs disappeared when *I*~Na~, *I*~f~, and *I*~K1~ were all switched to the LE values, although AP was inducible by external stimuli. The MDP gradually shifted to the negative direction and the overshoot potential became larger as *I*~Na~ increased, followed by the disappearance of *I*~f~ (Table [4](#Tab4){ref-type="table"}). The BCL was originally 510 ms in the EE model, and the increase in *I*~Na~ shortened the BCL to 340 ms, which was the shortest among the 3 representative spontaneous APs.Fig. 6Simulated tracings of action potentials (APs), half sarcomere length (hSL), and accompanying ionic currents. The relative densities of *I* ~Na~, *I* ~f~, and *I* ~K1~ were sequentially switched to the LE values from the EE model. An external stimulus (−38 pA/pF) was applied to the model in which the relative densities of *I* ~Na~, *I* ~f~, and *I* ~K1~ were all set to the LE valuesTable 4Characteristics of representative spontaneous action potentialsCombinationsMDP (mV)Overshoot potential (mV)BCL (ms)DSD (ms)EE model−81.402.14510300*I* ~Na~ set to LE−82.4837.16340130*I* ~Na~ and *I* ~f~ set to LE−84.7353.66780520Values were determined from the last spontaneous AP of the 600-s simulationMDP maximum diastolic potential, *BCL* basic cycle length, *DSD* diastolic slow depolarization

In the EE model, *I*~CaL~ was responsible for rapid depolarization of the membrane to overshoot the potential. As we sequentially switched *I*~Na~, *I*~f~, and *I*~K1~ to the LE values, the peak amplitude of *I*~CaL~ decreased from approximately −3.57 to −0.38 pA/pF and that of *I*~Na~ increased from 0 to −178.57 pA/pF, and *I*~Na~ became responsible for rapid depolarization rather than *I*~CaL~. Although we observed variations in the inward currents responsible for rapid depolarization, there were only slight differences in the outward currents (*I*~Kr~, *I*~K1~, and *I*~NaK~) responsible for membrane repolarization.

Discussion {#Sec14}
==========

We predicted membrane excitation patterns from among computer simulations of 512 combinations with the Kyoto model (Suppl. Fig. 1) and confirmed the simulated results with 2 other models, the TP model (Suppl. Figs. 2 and 3) and the LRd model (Suppl. Fig. 4). In all 3 models, \[Ca^2+^\]~i~ was increased to nonphysiological levels when *I*~K1~ was increased to the LE value before the disappearance of *I*~f~.

Burst-like membrane potentials in the Kyoto and TP models {#Sec15}
---------------------------------------------------------

Burst-like APs were observed in 8 combinations when the relative densities of *I*~st~ and *I*~CaL~ were set to the LE values and those of *I*~Na~, *I*~f~, *I*~K1~, and *I*~NaCa~ were set to the EE values (Suppl. Fig. 1); the burst-like APs disappeared as one of *I*~Na~, *I*~f~, *I*~K1~, and *I*~NaCa~ densities was shifted to LE value. In the TP model, on the other hand, burst-like APs were observed in 8 combinations in which the relative densities of *I*~f~ and *I*~NaCa~ were set to LE values and those of *I*~K1~ and *I*~st~ were set to EE values, and 2 combinations in which the relative densities of *I*~Na~, *I*~f~, *I*~CaL~, *I*~NaCa~, and *I*~st~ were set to LE values, and that of *I*~K1~ was set to the EE value (Suppl. Fig. 3). Although the combinations in which burst-like APs were observed were completely different between the Kyoto and the TP model, we observed similar dynamic changes in \[Na^+^\]~i~, which was increased during the repetitive bursts and decreased during the quiescent state between bursts (Fig. [3](#Fig3){ref-type="fig"}). During repetitive bursts, the amplitude of *I*~NaK~ was increased, which then decreased the amplitude of APs. As the bursts were terminated, *I*~NaK~ gradually decreased and contributed to the gradual increase in membrane potential during the quiescent state between bursts.

A similar pattern exhibiting burst-like APs (Fig. [3](#Fig3){ref-type="fig"}) has been reported in the pulmonary vein of rodents \[[@CR27]\], and such APs in the pulmonary vein are known to cause atrial fibrillation \[[@CR28], [@CR29]\]. In our simulations, burst-like membrane potentials were observed when the relative densities of *I*~st~ and *I*~CaL~ were set to the LE values and those of *I*~Na~, *I*~f~, and *I*~K1~ were set to the EE values in the Kyoto model. *I*~st~ has been reported as ionic currents of Na^+^ and K^+^ and has been observed only in SAN cells \[[@CR30]\]; however, there is no evidence for the presence of *I*~st~ in embryonic ventricular cells. Furthermore, we have no evidence for burst-like activities in ventricular cells at any stages of development. Although it may be interesting to note that such burst-like APs were predicted using both the Kyoto and TP models, we cannot draw notable conclusions from the simulated results because the combinations in which burst-like APs were observed were completely different between the Kyoto and TP models.

*I*~f~ should disappear before the increase in *I*~K1~ to avoid high intracellular Na^+^ and Ca^2+^ concentrations {#Sec16}
------------------------------------------------------------------------------------------------------------------

We observed abnormally high \[Ca^2+^\]~i~ values in all three models when *I*~K1~ was increased before the disappearance of *I*~f~; the simulated \[Ca^2+^\]~i~ was 12.1 μM in the LRd model, 6.27 μM in the TP model, and oscillated between 0.46 and 3.9 μM in the Kyoto model, when the relative *I*~K1~ density was set to LE and all the other densities were set to EE. In the TP model, the outward K^+^ current was increased as the relative *I*~K1~ density was increased from 0.11 to 1.0. Although the increase in outward K^+^ current polarized the membrane, the net inward current through *I*~f~ depolarized the membrane; therefore, the RMP was approximately −62.9 mV when the relative *I*~K1~ density was increased before the disappearance of *I*~f~. The increase in the inward Na^+^ current through *I*~f~ subsequently increased \[Na^+^\]~i~, which decreased the amount of Ca^2+^ excluded from the cytoplasm through *I*~NaCa~. Therefore, the high \[Ca^2+^\]~i~ in our simulation is mostly attributed to the decrease in *I*~NaCa~.

In the Kyoto model, the high \[Ca^2+^\]~i~ further led to the activation of *I*~LCCa~ (Fig. [4](#Fig4){ref-type="fig"}c), a current whose open probability increases at high \[Ca^2+^\]~i~ and contributes to transient inward current \[[@CR31]\]. The activated *I*~LCCa~ contributed to the depolarization of the membrane to −35 mV; subsequent activation of *I*~K1~ resulted in repolarization of the membrane to −80 mV. Although abnormally high \[Na^+^\]~i~ and \[Ca^2+^\]~i~ were observed in all three models, the spontaneous oscillation with long AP duration between −35 and −80 mV was observed only in the Kyoto model; therefore, we should note that such spontaneous oscillations are unlikely to be observed in actual embryonic ventricular cells.

In addition to the abnormally high \[Ca^2+^\]~i~, we observed that the RMP was slightly positive to −80 mV when *I*~K1~ increased before the disappearance of *I*~f~; the average RMP of 8 combinations in the LRd model was −69.1 mV and that of 32 combinations in the TP model was −62.9 mV. Of the 128 combinations in which *I*~K1~ was increased before the disappearance of *I*~f~ in the Kyoto model, the membrane potentials did not oscillate but were quiescent at approximately −50 mV in 32 combinations for which the relative *I*~CaL~ density was set to the EE value and the relative *I*~st~ density was set to the LE value. None of the combinations could fire APs even after the application of a large external stimulus in all 3 models because the depolarized membrane caused voltage-dependent inactivation of *I*~Na~, whereas abnormally high \[Ca^2+^\]~i~ caused Ca^2+^-dependent inactivation of *I*~CaL~. Although there is no evidence to suggest that quiescent cells from 12-day fetal rat hearts failed to evoke APs upon the application of external stimulus, it is worth noting that the predicted RMPs were roughly consistent with the experimental observations obtained from 12-day fetal rat hearts; Nagashima et al. \[[@CR5]\] obtained both cells with spontaneous APs and quiescent cells from 12-day fetal hearts wherein the quiescent cells exhibited an RMP of −48.4 ± 1.8 mV, which is more positive than the RMP of quiescent cells from 18-day fetal hearts (−80.9 ± 1.8 mV) \[[@CR5]\].

Relative densities of *I*~Na~ and *I*~f~ determine the BCLs and the overshoot potentials of the regular spontaneous APs {#Sec17}
-----------------------------------------------------------------------------------------------------------------------

Of the 9 components shifted between EE and LE values, the developmental changes in *I*~Na~ and *I*~f~ had large contributions to the variation in the BCL and the overshoot potential in the Kyoto model; the BCLs of the 160 regular spontaneous APs ranged from 306 to 884 ms and the overshoot potential ranged from 1.6 to 54.3 mV. We further shifted the relative densities of *I*~Na~ and *I*~f~ independently by a 10 % increment from EE to LE values and showed that the increase in *I*~Na~ shortened the BCL, whereas the decrease in *I*~f~ prolonged the BCL; the BCLs of the regular spontaneous APs were prolonged up to 1,550 ms when the relative *I*~Na~ density was set to 0.07 (EE value) and the relative *I*~f~ density was set to 0.1 (90 % shift to LE).

EE ventricles have a large range of BCLs---337--542 ms in 9.5-dpc mice \[[@CR6], [@CR8], [@CR9]\] and 273--2,500 ms in 12.5-dpc rats \[[@CR5], [@CR7]\]---and the beating rhythms of embryonic ventricular cells are irregular as reported in both 11.5-dpc rat ventricular cells \[[@CR5]\] and 18.5-dpc mouse ventricular cells with spontaneous APs \[[@CR6]\]. Although our simulations could not reproduce the irregular spontaneous APs reported in both mouse and rat embryonic ventricular cells, our predicted BCLs were all within the range of BCLs reported in experimental studies. We also showed that the wide range of the BCLs reported in vitro can be described by shifting the relative densities of *I*~Na~ and *I*~f~.

*I*~Na~ becomes responsible for membrane depolarization as *I*~Na~, *I*~f~, and *I*~K1~ are sequentially switched from EE to LE levels {#Sec18}
--------------------------------------------------------------------------------------------------------------------------------------

The simulated results imply that the increase in the relative *I*~K1~ density before the disappearance of *I*~f~ results in high \[Ca^2+^\]~i~ in all three models. In addition, we showed that the relative densities of both *I*~Na~ and *I*~f~ determine the BCL and overshoot potential of the regular spontaneous APs. On the basis of all the observations, we illustrated the representative changes in APs in which *I*~Na~ was increased before the disappearance of *I*~f~, followed by an increase in *I*~K1~ (Fig. [6](#Fig6){ref-type="fig"}). Following the sequence with representative models, we observed that *I*~Na~ took over the role of *I*~CaL~, which was originally the current responsible for the depolarization of the membrane in the EE model. This change in the dependence of depolarization from the Ca^2+^ current to the Na^+^ current is consistent with experimental observations in rodent ventricular myocytes \[[@CR32]\] in which the MDP shifted to a negative direction, also consistent with our simulation (Table [4](#Tab4){ref-type="table"}). The changes in BCL were roughly consistent with experimental observations on rat embryonic hearts; the BCL of the proximal ventricle in the 11.5-dpc embryonic rat heart was shorter than that of the ventricle in the 10.5-dpc rat; however, the BCL was prolonged again in the 12.5-dpc rat \[[@CR7]\].

Our hypothesis that an increase in *I*~Na~ density and disappearance of *I*~f~ should be observed in the early stage of embryonic development is supported by experimental observations that the densities of *I*~Na~ and *I*~f~ change earlier than those of other components \[[@CR3], [@CR6]\], including *I*~CaL~, *I*~K1~, *I*~NaCa~, and SR-related components \[[@CR2], [@CR20]\]. We demonstrate here that switching all the components in the mathematical model enabled us to simulate all possible combinations and identify pivotal component switches to describe the reported characteristics of embryonic ventricular cells. Our simulation procedure, together with experimental observations in the literature, will likely be useful in identifying the sequential regulation of gene or protein expression during development, which is difficult to determine through experimental data alone.

Limitations {#Sec19}
-----------

Our study has several limitations. The densities of ionic components were obtained from various rodents, including rats, mice, rabbits, and guinea pigs (Tables [1](#Tab1){ref-type="table"}, [2](#Tab2){ref-type="table"}, [3](#Tab3){ref-type="table"}), and implemented in the Kyoto and LRd models, both of which represent guinea pig ventricular cells. Therefore, simulations with the Kyoto and LRd models represent electrical activities of rodent ventricular cells in general. Although the TP model represents human ventricular cells, we implemented the densities listed in Tables [1](#Tab1){ref-type="table"}, [2](#Tab2){ref-type="table"}, and [3](#Tab3){ref-type="table"} in order to confirm our simulations with the Kyoto model; thus, the simulation with the TP model is not intended to represent developmental changes in human embryonic ventricular cells. In addition, changes in mRNA subtypes of the genes encoding the *I*~f~ and *I*~K1~ currents were not considered; the gene responsible for *I*~f~ is known to switch from *HCN4* to *HCN2* \[[@CR6]\], whereas that for *I*~K1~ switches from *Kir 2.2* to *Kir 2.1* \[[@CR5]\] during embryonic development. Although the length of the sarcomere was adopted as an index for the force of contraction, we did not consider developmental changes in contractile proteins; α- and β-myosin heavy chains (MHC), for example, are coexpressed and equally abundant in early embryonic ventricular cells, but α-MHC becomes predominant in adult ventricular cells \[[@CR33]\]. Therefore, the simulated changes in the length of the sarcomere may not be quantitatively accurate.

Conclusions {#Sec20}
===========

The relative densities of ionic components in mathematical models were switched independently between the EE and LE stages to identify pivotal components to describe reported characteristics of embryonic ventricular cells; all simulations were conducted using three models, the Kyoto, TP, and LRd models. In all three models, our simulations suggested that the tenfold increase in *I*~K1~ before the disappearance of *I*~f~ results in abnormally high \[Ca^2+^\]~i~. The developmental changes in relative densities of *I*~Na~ and *I*~f~ had large contributions to the wide range of BCL values in the regular spontaneous APs. Of the remaining six components in the Kyoto model, increases in *I*~CaL~ and SR-related components were involved in the enhancement of cell contraction.

Electronic supplementary material
=================================

 {#Sec21}

###### 

**Supplemental Figure 1.** Simulated membrane potential and half sarcomere length (hSL) for the 512 combinationsThe relative densities the 9 components, Na^+^ current (*I* ~Na~), funny current (*I* ~f~), inward rectifier K^+^ current (*I* ~K1~), sustained inward current (*I* ~st~), L-type Ca^2+^ current (*I* ~CaL~), Na^+^/Ca^2+^ exchange current (*I* ~NaCa~), ATP-sensitive K^+^ current (*I* ~KATP~), transient outward current (*I* ~to~), and sarcoplasmic reticulum (SR)-related components were switched independently to yield 512 combinations. The membrane potentials and hSLs are represented along the horizontal and vertical axes, respectively; the ranges of the axes are indicated in graphs. An external stimulus (38 pA/pF) was applied to the combinations shown as*blue* hysteresis loops at a frequency of 2.5 Hz.*Black letters in white boxes* in the column and row headers indicate that the relative density of the current was set to the early embryonic (EE) value, and*white letters in black boxes* indicate that the relative density of the current was set to the late embryonic (LE) value. (PDF 5332 kb)

###### 

**Supplemental Figure 2.** Simulated dynamics of intracellular Ca^2+^ concentration in 128 combinations using the Ten Tusscher--Panfilov model. The relative densities of *I* ~Na~, *I* ~f~, *I* ~K1~, *I* ~CaL~, *I* ~st~, *I* ~NaCa~, and SR-related components were switched independently between the EE and LE values. The*red* traces are illustrated from the time the first maximum diastolic potential (MDP) appeared after 600-s simulations. An external stimulation (−52 pA/pF at 1.0 Hz) was applied to the combinations without spontaneous activities. The paced simulations were conducted for additional 600 s to produce APs, represented as*blue* traces, which are illustrated from 600 to 602 s after the additional 600-s paced simulations. (PDF 858 kb)

###### 

**Supplemental Figure 3.** Simulated membrane potentials for 128 combinations using the Ten Tusscher--Panfilov model. The relative densities of *I* ~Na~, *I* ~f~, *I* ~K1~, *I* ~CaL~, *I* ~st~, *I* ~NaCa~, and SR-related components were switched independently between the EE and LE values. The*red* traces are illustrated from the time the first maximum diastolic potential (MDP) appeared after 600-s simulations. An external stimulation (−52 pA/pF at 1.0 Hz) was applied to the combinations without spontaneous activities. The paced simulations were conducted for additional 600 s to produce APs, represented as*blue* traces, which are illustrated from 600 to 602 s after the additional 600-s paced simulations. (PDF 699 kb)

###### 

**Supplemental Figure 4.** Simulated membrane potentials for 32 combinations using the Luo--Rudy model. The relative densities of *I* ~Na~, *I* ~CaL~, *I* ~f~, *I* ~K1~, and SR-related components were switched independently between the EE and LE values. Spontaneous action potentials (APs), represented as*red* traces, were observed when the relative density of *I* ~K1~ was set to the EE value; the*red* traces are illustrated 1 s from the time the first maximum diastolic potential (MDP) appeared after 600-s simulations. Spontaneous APs disappeared when the relative density of *I* ~K1~ was set to the LE value, and external stimulus (−80 pA/pF at 2.5 Hz) was applied for additional 600 s to produce APs, represented as*blue* traces, which are illustrated from 600 to 600.8 s after the additional 600-s paced simulations. (PDF 471 kb)

This research was supported by funds from the Yamagata Prefectural Government and Tsuruoka City, Japan. We would like to thank the members of WGSP at the Institute for Advanced Bioscience, Keio University, for critical suggestions.

Conflict of interest {#d30e4002}
====================

The authors declare that they have no conflicts of interest.
