Mechanisms Underlying Spontaneous Action Potential Generation Induced by Catecholamine in Pulmonary Vein Cardiomyocytes: A Simulation Study

Abstract

Cardiomyocytes and myocardial sleeves dissociated from pulmonary veins (PVs) potentially generate ectopic automaticity in response to noradrenaline (NA), and thereby trigger atrial fibrillation. We developed a mathematical model of rat PV cardiomyocytes (PVC) based on experimental data that incorporates the microscopic framework of the local control theory of Ca²⁺ release from the sarcoplasmic reticulum (SR), which can generate rhythmic Ca²⁺ release (limit cycle revealed by the bifurcation analysis) when total Ca²⁺ within the cell increased. Ca²⁺ overload in SR increased resting Ca²⁺ efflux through the type II inositol 1,4,5-trisphosphate (IP₃) receptors (InsP₃R) as well as ryanodine receptors (RyRs), which finally triggered massive Ca²⁺ release through activation of RyRs via local Ca²⁺ accumulation in the vicinity of RyRs. The new PVC model exhibited a resting potential of −68 mV. Under NA effects, repetitive Ca²⁺ release from SR triggered spontaneous action potentials (APs) by evoking transient depolarizations (TDs) through Na⁺/Ca²⁺ exchanger (AP_TDs). Marked and variable latencies initiating AP_TDs could be explained by the time courses of the α1- and β1-adrenergic influence on the regulation of intracellular Ca²⁺ content and random occurrences of spontaneous TD activating the first AP_TD. Positive and negative feedback relations were clarified under AP_TD generation.

1. Introduction

Several types of atrial fibrillation may be attributed to the ectopic activity of myocardial cells in the sleeves of pulmonary vein cardiomyocytes (PVCs) under augmented sympathetic stimulation [1,2,3,4,5,6]. Supporting this hypothesis, electrophysiological and histochemical experiments of rat PVCs by Okamoto et al. demonstrated generations of spontaneous action potentials under the influence of noradrenaline (NA) in dissociated PVCs [7,8] as observed in the tissue preparations [9].

The arrhythmogenic influencing of the stimulation of the α1-adrenergic receptor (AR) has been extensively studied in atrial cardiac myocytes [10,11]. These studies have indicated the functional significance of the co-localization of type II inositol 1,4,5-trisphosphate (IP₃) receptors (InsP₃Rs) with type II ryanodine receptors (RyRs). The Ca²⁺ release through InsP₃Rs from the sarcoplasmic reticulum (SR) was suggested to directly evoke a larger Ca²⁺ release via RyRs to initiate Ca²⁺ sparks. If multiple Ca²⁺ sparks overlap others, the resultant Ca²⁺ transient can trigger the temporal depolarization (TD) of the resting membrane through the Na⁺/Ca²⁺ exchanger (NCX), thereby triggering action potentials (AP_TDs). Co-immunostaining of atrial myocytes with antibodies against type II RyRs and type II InsP₃Rs revealed that they were localized on the SR membrane facing the subsarcolemmal space. This coupling of InsP₃R with RyR has scarcely been observed in the free running network of SR, where RyRs are expressed at a lower density. Thus, the excitation–contraction (E–C) coupling in atrial myocytes may be activated if InsP₃R is activated by the α1-AR stimulation. Spontaneous TD might occur at the junction of the SR terminal cisterna with the cell surface membrane in atrial myocytes.

In contrast to atrial myocytes, the E–C coupling in ventricular myocytes mostly occurs at the dyadic junction all over the cell, but the density of the InsP₃R expression is much lower than that in the atrial myocytes [10]. The density of the distribution of RyRs in ventricular myocytes is much higher than in the atrial myocytes. Therefore, the probability of direct coupling of RyRs with InsP₃R might be expected to be much less in ventricular myocytes, where the enhancement of TD is barely evoked by the α1-AR stimulation.

The PVCs in the sleeves of pulmonary veins are comparable in size and shape to ventricular myocytes and have a regular sarcomere pattern. The localization of InsP₃Rs was demonstrated to be parallel to the sarcomere pattern at a high density by Okamoto et al. [7]. They also demonstrated that the noradrenaline (NA) stimulation evoked a train of AP_TD discharge and that the resting membrane potential showed sporadic micro-fluctuations even in the absence of adrenergic stimulation.

Although different cell types might be involved in the arrhythmogenic activity in the pulmonary vein in different species [3], we focused on a particular cell type, which has been well characterized in a rat pulmonary vein [7,8]. We examined mechanisms underlying the generation of a train of spontaneous AP_TDs in response to NA stimulation by developing a mathematical PVC model. Until now, mathematical modeling of PVC has only been completed in humans [12] and rabbits [13], both of which showed automaticity in the absence of NA stimulation. Conversely, neither cardiomyocytes nor tissue preparation in PVs exhibited automaticity in rats [7,9]. To understand the mechanisms underlying the automaticity induced by NA stimulation by involving mobilization of intracellular Ca²⁺ as suggested by Okamoto et al., we used a refined model of the Ca²⁺-induced Ca²⁺ release (CICR) mechanisms, which was developed based on the Hinch model of Ca²⁺ releasing unit CaRU [14,15] and was adjusted to the human ventricular cell (HuVEC) model [16]. The novel CICR model was adopted in our study in combination with the InsP₃R model. The desensitization of the response to continuous α1-AR stimulation [17] and β1-AR stimulation [18] were defined by empirical equations. The simulation results showed that the arrhythmogenic activity (spontaneous generation of repetitive AP_TD) of PVCs could be initiated when variable conditions are combined in addition to the close coupling of InsP₃R and RyR. The pivotal additional factor is the activation of the β1-AR receptor to increase the Ca²⁺ store within SR (Ca²⁺ overload) through enhancement of both L-type Ca²⁺ channel (LCC) and sarcoplasmic/endoplasmic reticulum calcium pump (SERCA). The whole cell simulation of the PVC model suggests that Ca²⁺ extrusion by NCX, the temporal Ca²⁺ accumulation near the Ca²⁺ releasing site within the cytosol, and a Ca²⁺ overload of the SR are also essential in evoking spontaneous AP_TDs in the PVCs.

2. Results

2.1. Intracellular Ca²⁺ Dynamics Revealed by Bifurcation Analyses

Experimental studies have established that spontaneous Ca²⁺ release from SR occurs, and the Ca²⁺ transient evokes TD through the activation of NCX if the cardiac myocyte is overloaded with Ca²⁺ [19,20,21,22]. These cytosolic Ca²⁺ oscillations are important because they are accompanied by spontaneous oscillations in current and membrane potential to produce action potentials [23]. To clarify the behavioral characteristics of the cytosolic [Ca²⁺] in our PVC model, we examined steady-state solutions in the minimal model of intracellular Ca²⁺ dynamics. The constituents of this minimal model are described in Section 4. In short, we removed all the current systems on the plasma membrane from the PVC model to focus on the oscillatory behavior of intracellular Ca²⁺. To this minimal model, we applied bifurcation analyses, which have been successfully used to identify equilibrium points and/or the stable limit cycle in the mathematical whole cell or reduced models [24,25,26].

Firstly, we fixed the total Ca²⁺ content (Ca_tot) within a cell at a control level (3.20 femtomole) and found that the system was stable at a Ca_tot at junction space (jnc; [Ca²⁺_tot]_jnc) of 0.143 mM. The stable point was on the red curve on the left in Figure 1A (upper panel). Then, the amount of total Ca²⁺ was used as a bifurcation parameter and increased from 3.20 to 7.20 femtomole. The result shown in Figure 1A revealed that the system diverged at Ca_tot = 4.180 femtomole. When Ca_tot increased over the diversion point in the range of 4.180 to 5.923 femtomole, stable limit cycles (green) appeared with unstable equilibrium points (black). These results indicate that cyclic Ca²⁺ release occurred when Ca_tot was in the range plotted in green. In the usual numerical integration simulation, the peak and bottom level of the Ca²⁺ transients agreed well with the upper and lower level of the limit cycle depicted in green. As the α1-AR stimulation is mediated by the InsP₃R, its open probability of the channel (pO_InsP3R) was used as a bifurcation parameter in Figure 1B. Ca_tot was fixed at 3.6 femtomole for this analysis. At this low Ca_tot level, a stable equilibrium resting [Ca_tot]_jnc = 0.178 mM was observed. When pO_InsP3R increased, stable limit cycles were evoked between pO_InsP3R 0.0047 and 0.067. Increasing pO_InsP3R induced stable limit cycles at a Ca_tot level at which the system was stable without InsP₃R activation, and vice versa, increasing Ca_tot induced stable limit cycles even when there was no InsP₃R activation in this system. As is evident from the lower panels in Figure 1, the frequency of the automaticity increased markedly with increasing Ca_tot or pO_InsP3R.

The connection between the limit cycle identified in the bifurcation analysis and the cycle of spontaneous [Ca²⁺] fluctuations in the bulk cytosol space (blk; [Ca²⁺]_blk) was examined by conducting numerical time integration of the same minimal model (Figure 1C). Around the threshold Ca_tot level for initiating the limit cycle in Figure 1A, the Ca_tot was increased with a step size of 0.02 femtomole. The minimal model was quiescent up to 4.16 femtomole, and the cyclic [Ca²⁺]_blk transient was firstly observed at 4.18 femtomole. The frequencies of the Ca²⁺ transient increased with increasing Ca_tot as shown in Figure 1A. These responses of the minimal model were completely reversible when Ca_tot decreased.

2.2. Absence of Automaticity Inherent in Plasma Membrane Ionic Channels in the PVC

The PVC model showed a resting potential (V_rest) at −68 mV, which was similar to that observed in experiments (>−75 mV). This less negative resting potential, compared to ventricular cells, is due to a much lower density of the inward rectifier potassium current (I_K1) distribution compared with ventricular myocytes. Okamoto et al. revealed a hyperpolarization-activated Cl^– current (I_Clh) at a voltage range more negative than the resting membrane potential (V_rest) [8]. This I_Clh potentially stabilizes the low V_rest by antagonizing any kind of hyperpolarizing current, such as the I_K1. I_Clh can provide the largest inward current during strong hyperpolarizing pulses, but at the V_rest level, the magnitude of I_Clh activation is the lowest in the present mathematical model (violet trace in Figure 2D lower panel). Although I_Clh was not calculated in most of the simulations, I_Clh potentially depolarized the V_rest by ~7 mV in the PVC model, as suggested in experiments by Okamoto et al. [8].

The membrane excitation was examined in a cell model after equilibration without application of electrical stimulation for ~100 s (Figure 2). The two records of action potential superimposed in Figure 2A were evoked by either a brief electrical stimulus (−5 pA/pF), or a longer and smaller pulse (190 ms in duration).

The AP duration was ~30 ms at −40 mV, which is in the same range of ventricular and atrial myocytes in rats [7,8,27], but is clearly shorter than AP in the SA node pacemaker cells in rats [28]. The maximum rate of AP increased was mediated by the transient component of sodium current (I_Na) and was ~98.6 V/s. Repolarization was mainly caused by the inactivation of I_Na, and was facilitated by the simultaneous activation of I_Kto and I_Kur (Figure 2C). The role of I_CaL in maintaining the plateau potential of cardiac AP was largely compromised by the much larger transient outward potassium current (I_Kto) and ultra-rapid outward potassium current (I_Kur). I_CaL triggered CICR from SR to evoke the peak amplitude of the Ca²⁺ transient of ~1 μM to induce the developed tension of the myofilaments (Figure 2B).

The record of membrane currents evoked by the long pulse are demonstrated in Figure 2C,D, which were used to examine membrane automaticity of the PVC. Depolarization to a less negative potential range than −60 mV evoked an AP, followed by brief damping oscillations, but failed to evoke repetitive APs. The V_m during the diastolic period is largely determined by the balance of major currents of outward-going I_NaK, the background currents of I_Kbg and I_Kr (Figure 2D, upper panel), and the inward-going I_Nabg and I_NCX (Figure 2D, lower panel). Even if larger depolarizing current pulses were applied, no repetitive APs were observed. We concluded that the inactivation of I_Na during the AP was not removed at the diastolic potential levels during the current pulse because the activation of I_Kto and I_Kur quickly ceased during the falling phase of the AP. The amplitude of I_Kr (depicted in green in Figure 2D) is smaller than the I_Kbg. This is in strong contrast to the pacemaker mechanism of SA node cells [28], where the I_CaL is relieved from inactivation through hyperpolarization caused by the activation of I_Kr by the preceding AP. The amplitude of I_Kr assumed in the present model (depicted in green in Figure 2D) was too small to generate an appreciable size of a tail current on repolarization, in agreement with the experimental study [8]. We concluded that no stable cycle of APs arises in the membrane ionic system of the PVC model.

2.3. Generation of AP_TDs and an Activation Threshold for Ca²⁺_tot

The bifurcation analysis (Figure 1) applied to the intracellular Ca²⁺ dynamics in the absence of membrane ion fluxes revealed that spontaneous and repetitive Ca²⁺ transients were initiated by increasing Ca_tot above a certain threshold level (Ca_tot,thr = 4.18 femtomole) in the absence of plasma membrane ion fluxes. In the simulation shown in Figure 3, the spontaneous activity was recorded in the presence of intact plasma membrane currents.

The intracellular Ca_tot was technically increased step by step by enlarging the scaling factor of I_Cabg. The spontaneous discharge of AP_TD appeared when Ca_tot was increased above 3.878 femtomole in the absence of AR stimulation, whereas the threshold Ca_tot decreased to 3.617 femtomole in the presence of 0.15 μM [ISO] and 0.15 μM [IP₃]. Although the cycle length was markedly different, in both cases, the time course of V_m was flat during diastole except for a small slow diastolic depolarization of a few mV observed over ~300 ms after the preceding AP (Figure 3A,C). The Ca²⁺ flux (J_Ca) from SR via RyR (passive leak component (J_{leak_SR}) and active release component (J_{rel_SR})), and via InsP₃R (J_InsP3R) recovered in parallel with the replenishment of [Ca²⁺]_SRrl within the following 400–500 ms because of the Ca²⁺ diffusion from SR uptake sites. The total continuous Ca²⁺ efflux (J_{rel_SR}) finally evoked the final full Ca²⁺ release as indicated by the rapid fall of [Ca²⁺]_SRrl through an increase in [Ca²⁺]_jnc.

To confirm the involvement of TD in triggering APs, the amplitude of TD was depressed by increasing the membrane I_K1 conductance temporarily by six-fold (Figure 3B) or by three-fold (Figure 3D) approximately 1–3 s before the start of the recorded segment in the figure. This intervention blocked the AP generation and revealed the presence of TD in response to the spontaneous Ca²⁺ release. The amplitude of TD was 10–13 mV and the V_m at their peak was ca. –61 mV, most probably more negative than the threshold potential of I_Na activation.

Note that the downward deflection of [Ca²⁺]_SRrl showed double peaks: One synchronized with the start of TD and the other with the foot of the AP (or the activation of I_CaL) every time when AP_TD was triggered, whereas the second peak disappeared when TD failed to trigger AP_TD.

The distribution of Ca_tot was 30% in SR uptake site (SRup), 24% in SR releasing site (SRrl), and 43% in blk mostly bound with the buffer, and the minor components were found in jnc (1.3%) and iz (2.3%) in control and remained within ±5% in various experimental conditions examined in the present study.

Alternatively, the AP_TDs could be evoked by applying NA at 0.15 μM at the control sfI_Cabg = 1.68. Figure 3C,D show AP_TD as in Figure 3A,B, and demonstrate that the interval between the two successive AP_TDs is much reduced (interval = ~500 ms), and the preceding TD evoked by the spontaneous Ca²⁺ release is obvious before the rising phase of the AP. Again, no slow diastolic depolarization was observed, and a rapid decay in J_Ca via a cluster of RyRs (couplon) clearly precedes the onset of TD, and the second notch of J_Ca was caused by the opening of LCC during the rising phase of the AP. The TD from the rising phase of the AP was isolated by increasing sfI_K1 by three-fold (Figure 3D). The decrease in the TD amplitude in Figure 3D compared to that in Figure 3B was due to the enlarged I_NaK through the accumulation of [Na⁺] during the AP_TD burst.

The threshold Ca_tot level for initiating the AP_TD was examined in four combinations of two [ISO] of 0 and 0.2 μM, with two [IP₃] of 0.015 and 0.15 μM. The threshold Ca_tot at the 0.15 μM [IP₃] was higher (3.86 and 3.83 femtomole) than the value of (3.65 and 3.59 femtomole) obtained at 0.015 μM [IP₃]. For a reference level of Ca_tot, a representative 3.86 and 3.59 femtomole levels will be indicated in graphs of plotting time courses of Ca_tot in the following section.

2.4. Separation of α1- and β1-AR Influences on Ca_tot under Resting Condition

The findings in both Figure 1 and Figure 3 indicate that the increase in Ca_tot is the key factor for the initiation of the AP_TD in the present PVC model. Therefore, activation of the individual target components LCC, SERCA, and the Na⁺/K⁺ pump via β1-AR stimulation were examined by modifying Ca_tot at a saturating concentration of ISO (0.2 μM). In Figure 4A, the increase in Ca_tot evoked by the full member activation of β1-AR stimulation (LCC, SERCA, and Na⁺/K⁺ pump) was compared with different combinations of target activation. The increase in Ca_tot evoked by activating LCC plus SERCA was slightly larger than the control response. The separated Na⁺/K⁺ pump activation (the last trace in Figure 4) showed that the initial slight increase in Ca_tot was soon converted to a decrease during the ISO application due to a gradual decrease in Ca_tot through the NCX-mediated Ca²⁺ extrusion driven by the accelerated Na⁺/K⁺ pump. The influence of LCC activation is synergistic with the SERCA activation because these two molecules work together in transporting Ca²⁺ from extracellular space to the SR space. Ca_tot was markedly decreased by increasing the cytosolic [IP₃] (Figure 4B), which was applied to represent the activation of α1-AR. The application of two concentrations of IP₃ caused a dose-dependent decrease in Ca_tot as shown by removing the desensitization of the receptor at [IP₃] of 0.075, and 0.15 μM. If the desensitization (DS) was intact (0.15 μM + DS), the initial decreasing phase of Ca_tot was interrupted by the time-dependent desensitization.

2.5. Marked Latency before the Onset of Repetitive AP_TD Generation after AR Stimulation

In both tissue [9] and isolated rat cell preparations [7], a train of spontaneous AP was evoked after an extraordinary long latency (several minutes) after the application of NA. A plausible mechanism was already demonstrated in Figure 4B. The transient decrease in Ca_tot induced by activating InsP₃R (α1-AR) might counteract the β1-AR activation, which promotes the AP_TD burst by increasing Ca_tot through I_CaL and SERCA activation (Figure 4A). This was the case, as shown in Figure 5(A2), where the Ca_tot rapidly decreased (by ~7%) after the start of AR stimulation. In this simulation, the control Ca_tot level was set slightly lower than the threshold of initiating AP_TD (green vertical line) in the absence of AR stimulation. Factors involved in this response are shown in Figure 5(A4). The open probability (pO_InsP3R) of InsP₃R (green) quickly increased from 0.00027 to a peak of 0.00385 and then decayed to 0.00039 due to the desensitization of the α1-AR receptor activation. The Ca_tot followed a time course similar to pO_InsP3R (Figure 5(A2)). The slightly higher desensitized pO_InsP3R than the control was due to a rise in Ca_tot. Note, the InsP₃R can be partially activated by Ca²⁺. Figure 5(A3) indicates that the [Ca²⁺]_SRup and [Ca²⁺]_SRrl, the major store site of Ca_tot, followed the similar time course as Ca_tot.

The time courses of the β1-AR stimulation of LCC (af_CaL, red trace in Figure 5(A4)), SERCA (af_SERCA, black), and Na⁺/K⁺ pump (af_NaK, blue) are biphasic because of delayed fractional (30%) desensitization of β1-AR. The discharge of AP_TD evoked a full CICR, and thereby [Ca²⁺]_SRrl decreased to a minimum level (Figure 5(B3)). The [Ca²⁺]_SRup increased approximately parallel to the rise in [Ca²⁺]_blk, and showed a transient increase due to the uptake of Ca²⁺ via SERCA at every Ca²⁺ transient. The wash out of the influence of AR stimulant took an exponential time course and lasted for ~2 min.

2.6. Involvement of the Spontaneous Membrane Fluctuations in Determining the Latency

In the control run shown in Figure 5A, the PVC remained quiescent because the Ca_tot level was still lower than the reference level (green line) of the spontaneous AP_TD discharge even after the desensitization of α1-AR. The AP burst was initiated only when the Ca_tot was further increased by augmenting the β1-AR stimulation. In such case, a train of AP_TDs started at a given delay, different from the experimental finding of large variation in the 0–10 min latency [7]. For example, if scfI_Cabg was increased from 1.68 to 2.2, Ca_tot was still lower than the reference level, but the application of a saturating concentration of ISO (0.2 μM) evoked a train of AP_TD after a latency of 181.94 ± 0.5 s (mean ± SD, n = 13). This is in strong contrast to the marked experimental variations in the latency. We suggest variable time courses of gradual accumulation of Ca_tot after the AR stimulation. In the simulation (Figure 5(A2)), however, the accumulation of Ca_tot reached a peak within two minutes. We examined an alternative possibility of the miniature V_m fluctuations, most probably caused by sporadic CICR. The sporadic CICR was imitated using a random function (see Section 4) and the simulation results demonstrated in Figure 5B,C were obtained using the same protocol as in Figure 5A at the standard sfI_Cabg = 1.68 (Ca_tot = 3.765 femtomole). The frequency of sporadic CICR was approximately adjusted to the experimental records using a conventional RND function (see Section 4). The timing of random stimulation is shown by the vertical bars below the V_m record. Figure 5B,C provide typical examples of both success and failure of triggering repetitive AP_TDs. The size of the sporadic CICR was variable at individual CICR events, as evidenced by the miniature fluctuations in [Ca²⁺]_SRrl in the absence of continuous AP_TD generation (Figure 5C). Thus, the sporadic CICR mostly failed to evoke TD in visible amplitude to evoke AP_TD, and the AP_TD was triggered randomly. In the case of Figure 5C, two sporadic AP_TDs were evoked toward to the end of the AR stimulation, but failed to start the repetitive AP_TDs. In contrast, enough TD amplitude was evoked before the desensitization of β1-AR to trigger the AP_TD burst at a higher probability in Figure 5B.

The repetitive generation of AP_TDs was progressively accelerated as indicated by the plot of discharge frequency (Figure 5(B4), green curve). This is consistent with the bifurcation analysis results (Figure 1) that the frequency of the limit cycle of events increased with increasing Ca_tot in the Ca²⁺ dynamics when separated from the membrane function. In Figure 5B, the discharge of AP_TD increased Ca_tot through the activation of membrane I_CaL, and thereby positive feedback is involved in the rising phase of the spontaneous rate. This positive feedback was counteracted by the progressive increase of outward-going I_NaK through the accumulation of [Na⁺]_cyt (Figure 5(B4)), which was induced by the Na⁺/Ca²⁺ exchange via NCX. If the accumulation of [Na⁺]_cyt was augmented by increasing the β1-AR simulation, the increase in outward I_NaK current blocked triggering AP_TD by decreasing the amplitude of TD (not shown). In this case, the repetitive AP_TD generation was interrupted, but resumed when the accumulation of [Na⁺]_cyt was relieved by the extrusion of Ca²⁺ via NCX. Thus, the AP_TD burst occurred intermittently (not shown). Note, [Na⁺]_cyt was decreased via augmentation of the Na⁺/K⁺ pump by the β1-AR stimulation when the repetitive AP_TDs were absent (Figure 5A,C).

2.7. Latency Histogram of the AP_TD Burst

The latency was measured by constructing a latency histogram from several data sets of 1000 trials of the AR stimulation, as shown in Figure 6, where each data set was obtained at the control Ca_tot (scfI_Cab = 1.68) but at different amplitudes of I_CaL. In the control run (I_CaL × 1.0, Figure 6A), the burst was successfully evoked in 996 trials, and the histogram showed a peak at bin No. 6 of 2.5–3.0 min. The success rate evoking a repetitive AP_TDs was sensitive to the amplitude of I_CaL (Figure 6B–D), and the success rate drastically decreased to 33.4% by reducing the amplitude of I_CaL by 15% (Figure 6D). This success rate is near to that (26.7%) obtained in the experimental study [7]. The size of Ca²⁺ influx due to the activation of I_CaL and the uptake of Ca²⁺ into SR by SERCA play a key role in determining the time-dependent increase in Ca_tot. Note, the [Ca²⁺]_nd is largely dependent on [Ca²⁺]_SRrl when a continuous Ca²⁺ leak is generated by both InsP₃Rs and leak conductance of couplons.

3. Discussion

3.1. Brief Summary

The PVC model developed in the present study reconstructed well the electrical activity recorded in the PVCs dissociated from rat pulmonary veins by Okamoto et al. [7,8]. The model structure of CICR was refined in the human ventricular cell model [16], and was adopted in the presented PVC model including the Ca²⁺ spaces in the cytosol (Figure 7). The initiation of spontaneous APs in the model was attributed to the CICR (Figure 3), which was enhanced by the NA application, whereas the intrinsic membrane ionic mechanisms, described in cardiac pacemaker cells are not sufficiently developed to trigger spontaneous APs (Figure 2). After all membrane ionic fluxes were removed, the bifurcation analyses disclosed stable limit cycles over a certain range of Ca_tot within the cell (Figure 1). The threshold level for the initiation of the repetitive AP_TD generation was slightly lower in the full model than in the model of cytosolic Ca²⁺ dynamics (Figure 5), most probably because the membrane components, such as NCX and LCC, which enhanced the positive feedback mechanisms of the spontaneous rhythm, are not involved in the minimal model. We demonstrated that the β1-AR stimulation increased Ca_tot, whereas the α1-AR activation temporary decreased Ca_tot (Figure 4). The marked latency of several minutes for the start of repetitive AP_TD generation after the onset of stimulation could be explained by assuming the desensitization of α1-AR (Figure 5). Since the experimental findings on both the intracellular Ca²⁺ dynamics and signal transduction in rat PVC cells are still limited, the model proposed in this study may be revised in future studies. The presented simulation model, however, should provide a useful working hypothesis for conducting new experiments using dissociated cells.

3.2. Co-Localization of InsP₃R with RyRs in the Sub-Sarcolemmal Space Supporting Spontaneous CICR

The experimental studies in rat atrial myocytes demonstrated that a long lasting Ca²⁺ release from SR is induced by IP₃, and through the co-localization of InsP₃R with RyRs, the transient Ca²⁺ release is evoked [10]. Essentially the same mechanism was demonstrated when the InsP₃R was activated by endothelin-1 (ET-1) [11]. The Ca²⁺ extrusion through NCX is augmented by this transient increase in [Ca²⁺] and thereby a temporal transient membrane depolarization is evoked due to the electrogenic stoichiomety of Na⁺/Ca²⁺ exchange [19,20,21,22]. Mackenzie et al. directly recorded Ca²⁺ sparks as well as premature APs in the presence of ET-1 at a high probability of 75% in trials in rat atrial myocytes under α1-AR stimulation [11]. They demonstrated that a TD was evoked by an overlap of multiple Ca²⁺ sparks. The involvement of InsP₃R was proven by recording larger amplitudes of premature AP when a membrane-permeable analogue of IP₃ was applied in the extracellular medium [10]. The results of this simulation study using a new mathematical model of PVC in the present study agree with the basal Ca²⁺ mechanisms suggested by the previous studies. In the present study, a functional coupling of InsP₃R with RyRs was achieved by connecting the Ca²⁺ flux through InsP₃R to the hypothetical jnc, which directly encircles multiple CaRUs based on nanodomain (nd), consists of LCC and the couplon. Note that in the absense of LCC activation at the resting membrane potential, the spontaneous activation of a couplon is only controlled by the [Ca²⁺] within this jnc in the Hinch formalism of CaRU [14,15,16].

3.3. Peculiarity of AP_TD Generated in PVCs Compared to Aatrial and Ventricular Myocytes

The spontaneous TD and the AP_TD have been observed by activating the α1-AR more frequently in the atrial than in ventricular myocytes [10]. Atrial myocytes express InsP₃R at a much higher density than in the ventricular myocytes. InsP₃Rs are co-localized with the couplons (RyRs) in the junctional SR, but are absent in the network SR (non-junctional SR) in the atrial myocytes. The T-tubules are densely distributed at each interval of the sarcomere pattern, and thereby the sum of the T-tubule membrane is almost comparable to the surface membrane in ventricular myocytes (~78%) [29]. Thus, the junctional SR is distributed throughout the depth of the myocytes. However, the ratio of InsP₃R-coupled couplons to non-coupled couplons should be much lower in the ventricular myocytes. Thus, the spontaneous CICR may occur at a much lower frequency in the ventricular myocytes.

The PVCs show T-tubules and the InsP₃R distributed along the sarcomere pattern as in the ventricular myocytes [7]. Thus, it is expected that the density of couplons co-localized with InsP₃R should be much higher in PVCs than in atrial myocytes. This well agrees with PVs showing high frequencies of AP_TDs in both tissue preparations [2] as well as in dissociated cardiomyocytes [7], differently from the sporadic discharge of AP in the atrial myocytes during the activation of InsP₃R [10,11]. Note, the high frequency of AP_TD generation is attributable to the frequency of the intracellular Ca²⁺ oscillations as indicated by the bifurcation analysis (Figure 1). The low resting membrane potential (−70 mV) with low membrane conductance, most probably due to the low density of the I_K1 distribution, may further facilitate the discharge of repetitive AP_TDs. In conclusion, The AP_TDs may be evoked by essentially the same mechanism during AR stimulation as in the other working myocardial cells. The mechanisms may be highly enhanced in the PVCs to generate spontaneous repetitive AP_TDs in rat.

3.4. Ca²⁺ Overload and β1-AR Stimulation Evoking Repetitive AP_TD Generation in PVCs

Okamoto et al. observed that spontaneous TDs or miniature potential fluctuations occurred at the resting potential even before the application of NA [7]. These potential fluctuations are barely observed in the working myocytes, such as atrial or ventricular myocytes, but was induced with a high success rate by depressing the Na⁺/K⁺ pump by applying digitalis. This response is explained by the “Ca overload”, which was expected to occur secondarily to the increase in [Na⁺]_i [19,22]. The Ca²⁺ overload is also caused by the repetitive stimulation of ventricular myocytes in the presence of α1-NA stimulation, and the spontaneous Ca²⁺-release from SR evokes the TD [30]. We introduced a hypothetical I_Cabg for the purpose of varying the Ca_tot in PVCs. In the present study, the threshold of the AP_TD discharge was determined by slowly increasing Ca_tot in a stepwise manner to allow for a quasi-steady state distribution of Ca²⁺ within the cell in each step change. The threshold Ca_tot levels of 3.88, 3.87, 3.66, and 3.61 femtomole were determined in the control, α1-AR, β1-AR, and α1- plus β1-AR stimulation, respectively, after removing the desensitization of both receptors. This result is in good agreement with the finding in Figure 1 that the limit cycle of Ca²⁺ release was evoked simply by increasing the Ca_tot in the absence of the membrane ionic mechanisms.

3.5. The Latency and Frequency of AP_TDs Under NA Effects

The simulation results in Figure 5 strongly suggest that the latency before the initiation of the AP_TD generation might be determined by the desensitization of α1-AR in PVCs. In general, G_q-coupled receptors, such as α1-AR, ET-1 receptor, and angiotensin receptors, show desensitization after the binding of ligands [31,32,33]. The time course of desensitization has been demonstrated for the ET-1 and angiotensin receptors [34,35]. Cooling et al. [17] successfully developed a mathematical model of the signal transduction evoked by the G_q-coupled receptors according to the experimental measurements. They suggested that the time course of desensitization is largely determined at the level of receptor molecules. However, they did not discuss the desensitization time course of α1-AR. We failed to find any experimental measurements of the α1-AR desensitization time course in cardiac myocytes. We only observed several references in which indirect findings were described without interpretation by the authors. Zhang et al. recorded the I_CaL response to phentolamine application in rat ventricular myocytes, and found that the I_CaL amplitude temporarily decreased over several minutes before the gradual increase in I_CaL [36]. They suggested that the decrease might be due to Ca²⁺-mediated inactivation induced by a temporal release of SR Ca²⁺ through InsP₃R. Terzic et al. observed that the cell shortening evoked by the electrical stimulation of rat ventricular myocytes temporarily decreased in the same time course as I_CaL [37]. They suggested that this decrease in shortening might be due to the depletion of Ca²⁺ in SR. These findings strongly suggest the desensitization of α1-AR during phentolamine application. Although these responses were conducted at different ambient temperatures, the time courses were rather similar to the desensitization time course of α1-AR assumed in the present study.

In general, the subthreshold TDs appeared randomly as described in atrial myocytes [10,11]. The variation in the latent period before the initiation of the AP_TD generation is in agreement with the random occurrence of subthreshold TDs accompanied by transient elevations in intracellular Ca²⁺ [7]. To simulate this sporadic triggering of AP_TD, the regularity of spontaneous generation of AP_TD burst was avoided by lowering the Ca_tot to below the threshold level in the presence of β1-AR activation. Then, the train of AP_TD was triggered by the hypothetical random CICR. The experimental histogram obtained by Okamoto et al. [7] showed events of shorter latency than the time span of desensitization (~1.5 min). These shorter latency events might be simulated by increasing the conditioning Ca_tot, or by decreasing the conductance of whole-cell InsP₃R. In experiments, these parameters might be largely variable between individual PVCs, so that the latency histogram showed a wider distribution in the experimental study than in the present simulation study, in which the histogram was obtained using one cell model with a set of initial values of the variables.

3.6. Coupling Several Layers of Physiological Mechanisms to Regulate Ca_tot in PVCs

A regular time course of [Ca²⁺]_SRrl change (core cycle) is evoked when the cell is stimulated by an isolated electrical stimulus; the CICR transiently depletes [Ca²⁺]_SRrl at the terminal cysterna through excitation-contraction (EC) coupling. Then, the [Ca²⁺]_SRrl gradually recovers through the Ca²⁺ diffusion from the Ca²⁺ uptake site of SR. If Ca²⁺ leaks through RyRs and InsP₃R are added to this core cycle, any Ca²⁺ accumulation in jnc potentially triggers the next event of Ca²⁺ release via CaRU. The bifurcation analysis (Figure 1) indicated that a stable limit cycle of this Ca²⁺ release is evoked by increasing the Ca_tot beyond a certain level in the absence of Ca²⁺ flux across the cell membrane.

If the membrane ionic fluxes are coupled with the core cycle of Ca²⁺ oscillation, positive and negative feedback mechanisms are added. As a negative feedback, the transient Ca²⁺ release via couplons primarily accelerates extrusion of Ca²⁺ to the extracellular space via NCX. If an enlarged TD, evoked by the electrogenic Na⁺/Ca²⁺ exchange, triggers an AP_TD, the extra Ca²⁺ influx through LCC activation can increase Ca_tot, which may accelerate the core cycle of Ca²⁺ fluctuation and the spontaneous discharge of AP_TD to start the positive feedback cycle to further augment the Ca_tot.

The extrinsic regulation Ca_tot via β1-AR and α1-AR activation may further potentiate finer but complicated tuning of the Ca²⁺ dynamics. Basically, the simultaneous activation of both LCC and SERCA efficiently accumulates Ca_tot through Ca²⁺ uptake into SR. The activation of InsP₃R through α1-AR activation strongly decreases Ca_tot in combination with the Ca²⁺ extrusion by NCX. This simple antagonistic relationship between α1- and β1-ARs can be reversed to a synergistic relationship through the mechanism where by the frequency of the repetitive AP_TD generation is markedly increased by the faster recovery of [Ca²⁺]_SRrl by facilitating positive feedback.

So far, the influences of AR stimulation revealed in the present study were mostly suggested under the maximized strength of both α1- and β1-AR stimulation. Therefore, the mechanisms revealed by such pathophysiological simulations might not be applicable to the homeostatic regulation of Ca_tot under physiological conditions, where a much finer neural regulation may be expected. For example, the neural activity of peripheral sympathetic fiber always changes dynamically by repeating short bursts of APs in synchrony with fluctuations in arterial pressure [38]. Under such moderate stimulation, a vast desensitization, as observed in Figure 5 or in real experiments during the continuous stimulation for ~10 min, might not be expected. Most probably, the α1-AR stimulation may simply stabilize the positive inotropic influence through the β1-AR stimulation. The enhancement of Na⁺/K⁺ pump via β1-AR stimulation as well as the InsP₃R activation via α1-AR should compromise the increase in Ca_tot through an enhancement of NCX. The conduction of the AP from the atrial tissue may induce stable positive inotropy of PVCs through the β1-AR stimulation.

However, the simulation of the experimental findings successfully demonstrated cellular mechanisms in the present study. Thus, the mathematical model proposed is feasible for further examination of physiological mechanisms of the PVCs or might be applied to other cell types.

3.7. Limitations

The functional mechanisms of individual components of the presented PVC model had not yet been clarified by conducting experiments in PVCs. To overcome this issue, the components of the model were developed based on the general reaction scheme, which has been established in a variety of cardiac cell types. This situation of limited experimental data is similar to that of developing human cardiac myocytes. However, mathematical models of human cardiac myocytes have been useful in providing new working hypotheses. The PVC model provided an assumption that the relatively long latency of several minutes should be proved by directly examining the desensitization of α1-AR. Variations in Ca_tot should be estimated or considered when the time course of NA effect is examined.

Experiments in both tissue and dissociated myocytes demonstrated that the AP_TD burst was terminated with a long latency of many minutes after selectively blocking the α1-AR [2,7,9]. However, our simulation study failed to reconstruct this finding because the α1-AR was already desensitized to the control level in the presented simulation before the application of the α1-AR blocker. To address this issue, the involvement of other factors should be considered, such as the activation of I_CaL by the α1-AR pathway and/or a decrease in membrane K⁺ conductance, which takes a relatively long time course to develop. The nature of membrane K⁺ conductance (I_Kbg) depressed through α1-AR stimulation has not yet been clarified. Alternatively, the desensitization of α1-AR might partially occur and some fraction remains intact. If so, the pharmacologic blockade of α1-AR might decrease the frequency of repetitive AP_TD generation as indicated by the relationship between InsP₃R conductance and the frequency (Figure 1); thereby, the positive feedback cycle between the increase in spontaneous frequency of AP_TD through the accumulation of Ca_tot might be blocked, resulting in the cessation of spontaneous activity. We could simulate this mechanism under some narrow spectrum of the combination of modifying Ca_tot, degree of desensitization, or the decrease in the membrane K⁺ conductance. The complex mechanisms underlying the pluripotent nature of α1-AR stimulation remain to be clarified [36,37,39,40,41,42,43,44,45,46].

4. Materials and Methods

4.1. Intracellular Ca²⁺ Compartments and Distribution of Ionic Channels and Transporters in the PVC Model

Figure 7 shows a schematic representation of the model structure of a rat PVC. The T-tubule membrane provides the counterpart of the dyadic junction for the CICR. According to the Hinch model of the CaRU [14,15], individual CaRU is composed of one or a few number of LCC on the T-tubule membrane facing the couplons (clusters of RyRs) on the SR membrane. For activation and inactivation, both couplons and LCC refer to the same local Ca²⁺ concentration in the nanodomain (nd) cleft depicted in green. Each CaRU is separated from the others to support the local control of CICR, but a moderate cooperativity of multiple CaRUs is conserved by a local Ca²⁺ domain called jnc, which allows a temporal accumulation of released Ca²⁺ as firstly described in the HuVEC model [16,47]. Ca²⁺ accumulated in jnc gradually diffuses to iz and then to blk, in which myofilaments are located. Increased Ca²⁺ in each compartment is partly extruded by NCX. Note, only two representative CaRUs are demonstrated among many numbers of CaRUs in Figure 7. All CaRUs share a single common Ca²⁺ uptake site of the network SR spread in the blk for computational simplicity.

The major compartment of SR is SR_up equipped with SERCA, and provides its extension (terminal cysterna) to form SR_rl as the other counterpart of the dyadic junction. The SR_rl site contains the Ca²⁺ binding protein, calsequestrin, to increase the capacity of Ca²⁺ release. The Ca²⁺ released via InsP₃R diffuses into the jnc. The ion channels and transporters are exposed to different [Ca²⁺] after the Ca²⁺ release, and their distribution to each compartment is given in the Supplemental Materials, but not shown in Figure 7. During the Ca²⁺ release, [Ca²⁺]_SRrl is depleted to halt the Ca²⁺ release through the couplons.

NCX is distributed to jnc, iz, and blk, and their activation evokes TD. If the peak of transient depolarization crosses the activation voltage threshold of I_Na, AP_TD is triggered. A series of AP discharge occurs when Ca_tot within the cell is increased above a threshold level. In the present study, the extent of Ca²⁺ overload was adjusted simply by magnifying the hypothetical I_Cabg assumed on the surface membrane. Simulating various mechanisms of the Ca²⁺ overload was beyond the scope of the present study. The myofilament contraction model of Negroni and Lascano (2008) was adopted to calculate the free Ca²⁺ concentration in the blk.

4.2. Relationship Between Local [Ca²⁺]_nd and Spontaneous Ca²⁺ Release

In the Hinch model, the [Ca²⁺]_nd is determined in different combinations of LCC and couplons under the assumption of instantaneous Ca²⁺ distribution [14,15].

$${[C{a}^{2+}]}_{nd}=\frac{{[C{a}^{2+}]}_{jnc}+f_R\times {[C{a}^{2+}]}_{SRrl}+f_L\times \frac{\delta V·e^{-\delta V}}{1-e^{-\delta V}}\times {[C{a}^{2+}]}_o}{1+f_R+f_L\times \frac{\delta V}{1-e^{-\delta V}}},f_R=0.31,f_L=0.014.$$

(1)

The opening of LCC or a couplon is approximately represented by the parameters f_R and f_L defined by Equations (2) and (3), respectively. f_L and f_R are zero when LCC and couplons are closed, respectively.

$$f_R=\frac{G_R}{G_{diff}},$$

(2)

$$f_L=\frac{G_L}{G_{diff}},$$

(3)

where G_L and G_R are Ca²⁺ conductivity through the LCC and couplons, respectively; and G_diff represents the conductivity of Ca²⁺ diffusion across the hypothetical border between nd and jnc. In the case of couplons, its activation is assumed to occur in an all or none manner; thus, G_R represents the limiting conductance of the couplon, and f_R provides the proportional coefficient in respect to [Ca²⁺]_SRrl. In the case of LCC, a single LCC is assumed so that the conductance G_L represents the limiting conductance of LCC, but is an exponential function of the membrane potential V_m as described by Equation (1). Note, f_L and f_R also provide proportional coefficients for [Ca²⁺]_nd in respect to [Ca²⁺ ]_o in the extra-cellular solution and [Ca²⁺]_SRrl, respectively.

An increase in [Ca²⁺]_nd is the sole factor in determining the spontaneous Ca²⁺ release via couplon. Then, [Ca²⁺]_jnc is calculated from the total amount of Ca²⁺ within the space ([Ca_tot]_jnc), which is determined by the four Ca²⁺ fluxes as represented by Equation (4).

$$\frac{d{[C{a}_{tot}]}_{jnc}}{dt}=-\frac{J_{Ca,m}+J_{Ca\_couplon}+J_{Ca\_InsP3R}-J_{Ca\_jnciz}}{V_{jnc}},$$

(4)

where V_jnc is the volume of jnc, and J_Ca,m, J_Ca-couplon, J_Ca-InsP3R, and J_Ca,jnciz are membrane Ca²⁺ fluxes in the jnc space, Ca²⁺ release through couplon and InsP₃R, and Ca²⁺ diffusion between the jnc and iz, respectively. Under a Ca²⁺-overload condition, the J_{Ca_rel} through the basal openings of the couplon or activated InsP₃R increases in proportion to the level of [Ca²⁺]_SRrl above the threshold level. The amplitude of Ca²⁺flux is determined by Equations (5) and (6).

$$J_{Ca,RyR}=G_{RyR,bg}·({[C{a}^{2+}]}_{SRrl}-{[C{a}^{2+}]}_{jnc}),$$

(5)

$$J_{Ca,IP3R}=G_{IP3R,bg}·({[Ca]}_{SRrl}-{[Ca]}_{jnc}).$$

(6)

For Ca-buffers in jnc, iz, and blk, see Supplemental Materials.

4.3. Ion Channels and Transporters

All ion channels and transporters included in the present PVC model are listed with references in the Abbreviations section. All the ion channels and transporters are distributed in each Ca²⁺ compartment (

Table 1. Distribution of plasma membrane channels and transporters.

Channels and Transporters	jnc	iz	blk
Ca2+ insensitive channels (INa, IKr, IK1, IKto, IKur, IClh, IKbg, INabg)	-	0.1	0.9
IPMCA	-	0.1	0.9
INaK	-	0.1	0.9
ICaL(LCC)	0.75	0.15	0.1
INCX	0.03	0.25	0.72

).

The Cl⁻ channel having a large conductance was described during hyperpolarization in the rat PVC cells [8]. However, including I_Clh in the PVC model is difficult, since no experimental report on Cl⁻ transporters is available for the PVC cell. To obtain Cl⁻ homeostasis, we need to balance the passive Cl⁻ flux with some Cl⁻ transporters. At present, the PVC model does not include I_Clh. In preliminary simulations, I_Clh was included only to estimate its contribution to the membrane potential. We found that the contribution of I_Clh to the membrane potential is relatively small at V_m less negative than −70 mV.

The equations of these ion currents and the transporters are given in Supplemental Materials. Only several currents with immediate significance for the spontaneous APD_TD are described below.

In rat PVC cells, the amplitude of I_Na was larger when compared with I_Na in the atrial tissue, and the activation range was shifted to left in the current–voltage (I–V) relationship [6,48]. No obvious difference was reported in the density and kinetics between PVC and LA [3]. Okamoto et al. demonstrated that the inactivation of I_CaL is faster in the rat PVC cells compared to ventricular cells [7]. In the CaRU of the HuVEC model, the locus of Ca²⁺-mediated inactivation of LCC is nd and the inactivation rate (k_oc) is described as:

$$k_{oc}=\frac{{[C{a}^{2+}]}_x}{K_L}·\frac{f{V}_{m,act}}{T_L}{K}_L=0.0044mM,T_L=147.51,$$

(7)

where fV_m,act is a parameter for the V_m-dependent activation parameter and is used to describe the dependence of LCC inactivation on the opening of the voltage gate, and K_L is used to represent the [Ca²⁺]-dependence in nd or other [Ca²⁺]. We assumed that 75% of I_CaL was connected to jnc, and 15% to iz, and the rest to blk. A slightly smaller K_L (K_L = 0.00154 mM) was used for the at PVC cell model to increase the Ca²⁺ sensitivity.

Virtually no obvious amplitude of I_Kr tail current has been recorded on the jump from the positive potential to the holding potential in the voltage clamp experiment in rat PVCs [7]. However, in rat atrial cells [49,50], ventricular cells [50], and canine PVCs [3], the existence of I_Kr was revealed at a relatively small density. We included a minimum size of I_Kr in the PVC model to provide a repolarization reserve for technical reasons to avoid instability of V_m, which was observed around –20 to 40 mV after the inactivation of both I_Kto and I_Kur. Activation of a transient outward current was observed at the onset of the depolarizing pulse in the rat ventricular myocyte model by Pandit et al. [51]. The I_Kur model was adopted from the mouse ventricular myocyte model by Bondarenko et al. [52].

The I_Clh, found by Okamoto et al., showed slow voltage-dependent activation on hyperpolarization with a time course similar to that of the hyperpolarization-activated non-selective cation current in the SA node cell [8]. Thus, we applied the kinetic equation of I_h simplified by reducing the number of states to C₁, C₂, and O states after optimizing the rate constants. The I_K1 in PVCs seems much smaller in amplitude compared to the ventricular cells to allow the relatively low resting potential in PVC cells [7].

The histochemical examination of the NCX by Okamoto et al. [7] disclosed the localization of NCX near the T-tubules. This finding was represented by the 28% distribution of NCX near the T-tubule (2% in jnc and 25% in iz) and the rest in blk, respectively (Table 1).

4.4. Simulation of NA Stimulation of PVCs

The mathematical model of the β1-AR stimulation developed by Saucerman et al. [18] was adopted after minor modification [53,54]. The reaction cascade that generates cytosolic cyclic AMP (cAMP) and the activation of PKA evoked by β1-AR agonist ISO is described in the Supplemental Materials.

Doisne et al. never observed spontaneous activity in isolated pulmonary veins of rat under basal physiological conditions, but observed that the application of norepinephrine induced automatic electric activity [9], in agreement with the contractile activity observed by Maupoil et al. [2]. Okamoto et al. observed the spontaneous repetitive AP generation when α1- and β1-ARs were stimulated in dissociated PVCs, and also recorded major ionic current components of the PVCs [7]. In reconstructing the experimental data of Okamoto et al. [7,8], we selected SERCA, LCC, and the Na⁺/K⁺ pump as key targets of the β1-adrenergic regulation, and InsP₃R and a kind of background K⁺ conductance for the α1-AR regulation. Although RyRs are also the target of the β1-adrenergic regulation, we did not calculate the modification of RyRs for simplicity, partly because the activation of the couplon occurred in an all-or-none manner during the CICR. Here, the “catecholamine effect” denotes the sum of these effects. The influence of individual affecters were examined by varying the activity of the corresponding target molecule.

4.4.1. Implementation of β1-AR Effects

We assumed two populations of the target molecules for each of Na⁺/K⁺ pump, SERCA, and LCC; one is the phosphorylated fraction (Fr_PKA) by the catalytic subunit of PKA and the other is the dephosphorylated fraction by several kinds of phosphatase (PPs). In a normalized scheme, the reaction of phosphorylation is described by Equation (8):

$$(1-F{r}_{PKA})\begin{array}{c}\stackrel{\alpha }{\to }\\ \underset{\beta }{\leftarrow }\end{array}F{r}_{PKA}.$$

(8)

The time course of Fr_PKA after the onset of β1-AR stimulation was calculated using Equation (9).

$$\frac{dF{r}_{PKA}}{dt}=\alpha ·(1-F{r}_{PKA})-\beta ·F{r}_{PKA},$$

(9)

where a common time constant (τ) of phosphorylation was assumed for the three kinds of target protein. Since the activation time course takes several tens of seconds after the application of β1-AR stimulant, a τ of 40 s was assumed for convenience. The forward rate α in Equation (8) is dependent on the concentration of catalytic subunit of PKA (cat), a forward rate constant k_cat = 0.0625/mM/ms was determined by model adjustment, and the backward rate constant β was obtained by Equations (10) and (11).

$$\alpha =k_{cat}·cat,$$

(10)

$$\beta =\frac{1}{\tau }-\alpha .$$

(11)

When a basal phosphorylation (base) of target protein caused by kinases other than PKA is assumed, the sum of phosphorylated active fraction (af_(t)) of a target protein is given as a sum of two components denoted as base and delta.

$$af(t)=base+delta·F{r}_{PKA(t)}.$$

(12)

The proportion of base/(base + delta) for Na⁺/K⁺ pump, SERCA, and LCC were model adjusted to 25%, 10%, and 43.5%, respectively.

4.4.2. Na⁺/K⁺ Pump

The Na⁺/K⁺ pump model used in the human ventricular cell (HuVEC) model [55] was adopted, since it was developed by referring to a wide variety of electrophysiological findings to apply the detailed thermodynamic framework of the original Na⁺/K⁺ pump model [56]. The β1-adrenergic regulation of the Na⁺/K⁺ pump model was introduced by implementing the involvement of phospholemman (PLM) as described by Despa et al. [57]. According to their study, the activation of Na⁺/K⁺ pump by the phosphorylation of PLM was represented by a decrease to 75% the control in the apparent Na⁺-dissociation constant (Kd_,Nai) of the cytosolic binding site. The decrease was calculated by applying a scaling factor sf to the original K_d in the four state model of Na⁺/K⁺ pump.

$${\overline{Kd}}_{Nai}=s{f}_{NaK}·K{d}_{Nai}.$$

(13)

To satisfy the thermodynamic constraint of:

$$\frac{k_1^{+}{k}_2^{+}{k}_3^{+}{k}_4^{+}{K}_{d,K_i}^2{K}_{d.N{a}_e}^3}{k_1^{-}{k}_2^{-}{k}_3^{-}{k}_4^{-}{K}_{d,MgATP}{K}_{d,K_e}^2{K}_{d.N{a}_i}^3}=K_{~MgATP}·e^{F{V}_m/RT},$$

(14)

where the Kd,_Ke for the extracellular binding site was increased by sf^−3/2 to obtain the increase in the I_NaK [55,57], a sf_NaK = 0.72 was set for the Na⁺/K⁺ pump relieved from the inhibitory action of PLM.

We assumed two populations of the Na⁺/K⁺ pump, with phosphorylated PLM fractions (activated fraction, af) and non-phosphorylated fraction (1 – af). The whole cell I_NaK was determined as a sum of control I_NaK and activated ${\overline{I}}_{NaK}$,

$$I_{NaK}=(1-af)·I_{NaK}+af·{\overline{I}}_{NaK}.$$

(15)

4.4.3. SERCA

The biophysical model of SERCA that we developed by modifying the Tran et al. model [58] was used after adding the β1-adrenergic regulation. SERCA belongs to the same family of the P-type ion pump as the Na⁺/K⁺ ATPase; namely, the activity of SERCA was inhibited by the non-phosphorylated phospholamban (PLB), and SERCA was relieved from this inhibition through phosphorylation of PLB during the β1-AR stimulation. The regulation by AR was defined by decreasing the apparent K_d for the cytosolic Ca²⁺ to half the control (sf_SERCA = 0.5) when the PLB was phosphorylated. The thermodynamic constraint was satisfied by applying the same sf to the rate k_-1, which was used as an adjustable parameter in the original model. The time course of activation was defined by Equation (12) and the amplitude of the Ca²⁺ flux into SR was calculated by Equation (17). The same time constant τ = 1 s was used for the β1-AR regulation of Na⁺/K⁺.

$${\overline{Kd}}_{Cai}=s{f}_{SERCA}·K{d}_{Cqi}.$$

(16)

The whole cell J_SERCA was determined as a sum of control J_{SERCA_c} and activated J_{SERCA_a}.

$${\mathrm{J}}_{\mathrm{SERCA}}=(1-af)·{\mathrm{J}}_{\mathrm{SERCA}\_\mathrm{c}}+af·{\mathrm{J}}_{SERCA\_a}.$$

(17)

4.4.4. LCC

The amplitude of I_CaL increased during the β1-AR stimulation without a marked change in its time course. Thus, the amplitude of I_CaL is magnified when noradrenaline is applied to the model cell. We assumed that the fraction of non-phosphorylated population of LCC is 56.5% and shows virtually zero conductance in the absence of β-adrenergic stimulation. Thus, the whole cell conductance of I_CaL is proportional to af_(t) (Equation (12)).

4.4.5. α1-Adrenergic Signal, [IP₃]

We failed to find an α1-AR model to incorporate in the PVC model. In this study, the time-dependent activation and desensitization of the receptor was simply represented by changes in the active second messenger [IP_3a] using a concentration [IP₃]_(t), and an inactivation (i) parameter of the α1-AR at time t. To obtain [IP₃ ], equal to the control 0.015 mM after the desensitization, a lower limit of 0.1 was applied to the range of $i_{\infty }$:

$$[I{P}_{3a}]=i_{(t)}·{[I{P}_3]}_{(t)},$$

(18)

$$i_{(t+dt)}=i_{\infty }-(i_{\infty }-i_{(t)})·\exp (-\frac{dt}{{\tau }_i}){\tau }_i=30000(m{s}^{-1}),0.1\le i_{\infty }<1.$$

(19)

The [IP₃]_(t) at time t is also given by an exponential function:

$${[I{P}_3]}_{(t+dt)}={[I{P}_3]}_{\max }-({[I{P}_3]}_{\max }-{[I{P}_3]}_{(t)})·\exp (-\frac{dt}{{\tau }_i{}_{p3}}){\tau }_{ip3}=5000(m{s}^{-1}).$$

(20)

4.4.6. InsP₃R

In the cardiac myocytes, the type II homologous isomer is the major component of InsP₃R [59] and is expressed on the SR membrane near the T-tubules [7]. Thus, in the model, the InsP₃R was exposed to jnc, so that its Ca²⁺ release plays a critical role in determining the bias level for the CICR. The open probability pO_InsP3R was calculated using the InsP₃R-II model [60]. The influence of Ca²⁺ on the pO_InsP3R was much less than that of IP₃ over the range of [Ca²⁺] during the Ca²⁺ transient. The resting level of [IP₃] was set to 0.015 μM according to Cooling et al. [17] so that J_InsP3R was virtually zero. The [IP₃] was increased 10-fold (0.15 μM) during the stimulation of α1-AR. The Ca²⁺ flux via InsP₃R channel was calculated using Equation (21).

$$J_{InsP3R}=P_{IP3R}·p{O}_{InsP3R}·({[C{a}^{2+}]}_{SRrl}-{[C{a}^{2+}]}_{jnc}).$$

(21)

4.4.7. Background K⁺ Current, I_Kbg

Doisne et al. observed a slow continuous depolarization of ~+21 mV during a 15 min application of NA to the rat PV tissue, which was larger than the +12 mV depolarization in the atrial tissue [9]. Similar depolarization was observed in response to phenylephrine (a selective agonist for α1-AR) in the rat atrial tissue preparation by Jahnel et al. [61]. In our simulation, the AR-mediated decrease of the background I_Kbg was set at 20%–30% with a time constant of 120 s to obtain the experimental time course of depolarization.

$$I_{K_{bg}}=s{f}_{K_{bg}}·G_{K_{bg}}·(V_m-E_K),$$

(22)

$$s{f}_{K_{bg}(t+dt)}=s{f}_{K_{bg}(t)}+(s{f}_{K_{bg}}{}_{\infty }-s{f}_{K_{bg}(t)})·\frac{dt}{{\tau }_{Kb}}{\tau }_{Kb}=80000(m{s}^{-1}).$$

(23)

4.5. Simulation of the Random Events of CICR

If [Ca²⁺]_jnc fluctuates over a range slightly lower than the threshold of full CICR activation, the membrane potential may vary randomly. However, in the conventional numerical integration of d[Ca²⁺]/dt, the average time course of [Ca²⁺]_jnc only changes smoothly and continuously. In simulation, the random activation of a couplon could be technically evoked by introducing a random function in the present Hinch format of CICR. For practical purposes, the magnitude of [Ca²⁺]_jnc was multiplied two-fold at a probability of 0.0007 when the random function implemented in the VB system provided a value larger than (1 – 0.0007) at every step of numerical integration. This additional fluctuation randomly induced a full-blown CICR shown in Figure 5. Note, this modified [Ca²⁺]_jnc is only used to determine the opening rate of a couplon, but does not interfere with the mass conservation of Ca²⁺ dynamics.

4.6. Changes in Membrane Potential and Ion Concentrations

A C_m of 115.3 pF was assumed in the PVC model according to the measurements of Okamoto et al. [7]. The AP was triggered by applying a current pulse of −5 pA/pF for 3 ms. Changes in V_m as well as the intracellular ion concentrations were calculated.

$$I_{Ca,tot\_jnc}=I_{Ca,CaL\_jnc}+I_{Ca,NCX\_jnc},$$

(24)

$$I_{Ca,tot\_iz}=I_{Ca,CaL\_iz}+I_{Ca,bg\_iz}+I_{Ca,NCX\_iz}+I_{Ca,PMCA\_iz},$$

(25)

$$I_{Ca,tot\_blk}=I_{Ca,CaL\_blk}+I_{Ca,bg\_blk}+I_{Ca,NCX\_blk}+I_{Ca,PMCA\_blk},$$

(26)

$$I_{Ca,tot}=I_{Ca,tot\_jnc}+I_{Ca,tot\_iz}+I_{Ca,tot\_blk},$$

(27)

$$\begin{array}{l}{I}_{Na,tot}=I_{\mathrm{Na}}+I_{\mathrm{Na},\mathrm{CaL}\_\mathrm{jnc}}+I_{\mathrm{Na},\mathrm{CaL}\_\mathrm{iz}}+I_{\mathrm{Na},\mathrm{CaL}\_\mathrm{blk}}+I_{\mathrm{Na},\mathrm{bg}}\\ +I_{\mathrm{Na},\mathrm{NaK}}+I_{\mathrm{Na},\mathrm{NCX}\_\mathrm{jnc}}+I_{\mathrm{Na},\mathrm{NCX}\_\mathrm{iz}}+I_{\mathrm{Na},\mathrm{NCX}\_\mathrm{blk}}\end{array}$$

(28)

$$\begin{array}{l}{I}_{K,tot}=I_{\mathrm{k},\mathrm{K}1}+I_{\mathrm{K},\mathrm{Kr}}+I_{\mathrm{K},\mathrm{Kto}}+I_{\mathrm{K},\mathrm{ur}}+I_{\mathrm{K},\mathrm{bg}}+I_{\mathrm{K},\mathrm{Na}}\\ +I_{\mathrm{K},\mathrm{CaL}\_\mathrm{jnc}}+I_{\mathrm{K},\mathrm{CaL}\_\mathrm{iz}}+I_{\mathrm{K},\mathrm{CaL}\_\mathrm{blk}}+I_{\mathrm{K},\mathrm{NaK}}\end{array}$$

(29)

The concentration changes of intracellular ions were calculated by the numerical integration of the membrane fluxes across the surface membrane as well as the SR membrane.

$$\frac{dC{a}_{tot\_jnc}}{dt}=(-\frac{I_{Ca,tot\_jnc}·C_m}{2·F}+J_{Ca\_rel}-J_{Ca\_jnciz}+J_{IP3})/V_{vol\_jnc},$$

(30)

$$\frac{dC{a}_{tot\_iz}}{dt}=(-\frac{I_{Ca,tot\_iz}·C_m}{2·F}+J_{Ca\_jnciz}-J_{Ca\_izblk})/V_{vol\_iz},$$

(31)

$$\frac{dC{a}_{tot\_blk}}{dt}=(-\frac{I_{Ca,tot\_blk}·C_m}{2·F}-J_{Ca\_SERCA}+J_{Ca\_izblk})/V_{vol\_blk},$$

(32)

$$\frac{dC{a}_{SRup}}{dt}=\frac{J_{Ca\_SERCA}-J_{Ca\_transSR}}{V_{vol\_SRup}},$$

(33)

$$\frac{dC{a}_{tot\_SRrl}}{dt}=\frac{J_{Ca\_transSR}-J_{Ca\_rel}-J_{IP3}}{V_{vol\_SRrl}},$$

(34)

$$\frac{dNa}{dt}=\frac{I_{Na,tot·}{C}_m}{F·V_{vol\_cyt}},$$

(35)

$$\frac{dK}{dt}=\frac{(I_{K,tot·}+I_{apl})·C_m}{F·V_{vol\_cyt}},$$

(36)

$$\frac{d{V}_m}{dt}=-(I_{Ca,tot}+I_{Na,tot}+I_{K,tot}+I_{apl}).$$

(37)

4.7. Bifurcation Analysis on the Ca²⁺ Dynamics Within the Cell in the Absence of Membrane Ionic Fluxes

The model equations were all translated into an ordinary differential equation (ODE) file. Initial values were obtained using a numerical simulation program. Then equilibrium points were calculated using xppaut [62]. When Hopf bifurcation points were obtained, limit cycles were calculated using a conventional computational tool “Auto” to plot a bifurcation diagram. Ca_tot (attomole) was composed of various compartments as described by Equation (38).

$$C{a}_{tot}={[C{a}_{tot}]}_{jnc}·Vo{l}_{jnc}+{[C{a}_{tot}]}_{iz}·Vo{l}_{iz}+{[C{a}_{tot}]}_{blk}·Vo{l}_{blk}+{[C{a}_{tot}]}_{SRrl}·Vo{l}_{SRrl}+{[C{a}^{2+}]}_{SRup}·Vo{l}_{up}(\mathrm{attomole}).$$

(38)

In the bifurcation diagram and the main text, the unit of Ca_tot is converted to femtomole from attomole.