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

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 Ca2+ release from the sarcoplasmic reticulum (SR), which can generate rhythmic Ca2+ release (limit cycle revealed by the bifurcation analysis) when total Ca2+ within the cell increased. Ca2+ overload in SR increased resting Ca2+ efflux through the type II inositol 1,4,5-trisphosphate (IP3) receptors (InsP3R) as well as ryanodine receptors (RyRs), which finally triggered massive Ca2+ release through activation of RyRs via local Ca2+ accumulation in the vicinity of RyRs. The new PVC model exhibited a resting potential of −68 mV. Under NA effects, repetitive Ca2+ release from SR triggered spontaneous action potentials (APs) by evoking transient depolarizations (TDs) through Na+/Ca2+ exchanger (APTDs). Marked and variable latencies initiating APTDs could be explained by the time courses of the α1- and β1-adrenergic influence on the regulation of intracellular Ca2+ content and random occurrences of spontaneous TD activating the first APTD. Positive and negative feedback relations were clarified under APTD generation.


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 3 ) receptors (InsP 3 Rs) with type II ryanodine receptors (RyRs). The Ca 2+ release through InsP 3 Rs from the sarcoplasmic reticulum (SR) was suggested to directly evoke a larger Ca 2+ release via RyRs to initiate Ca 2+ sparks. If multiple Ca 2+ sparks overlap others, the resultant Ca 2+ transient can trigger the temporal depolarization (TD) the PVC model to focus on the oscillatory behavior of intracellular Ca 2+ . 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 2+ 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 2+ 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 2+ 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 2+ 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 2+ 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 3 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 3 R activation, and vice versa, increasing Ca tot induced stable limit cycles even when there was no InsP 3 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 . 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 2+ content (Catot) within a cell at a control level (3.20 femtomole) and found that the system was stable at a Catot at junction space (jnc; [Ca 2+ 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 2+ 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 Catot = 4.180 femtomole. When Catot 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 2+ release occurred when Catot was in the range plotted in green. In the usual numerical integration simulation, the peak and bottom level of the Ca 2+ 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 InsP3R, its open probability of the channel (pOInsP3R) was used as a bifurcation parameter in Figure 1B. Catot was fixed at 3.6 femtomole for this analysis. At this low Catot level, a stable equilibrium resting [Catot]jnc = 0.178 mM was observed. When pOInsP3R increased, stable limit cycles were evoked between pOInsP3R 0.0047 and 0.067. Increasing pOInsP3R induced stable limit cycles at a Catot level at which the system was stable without InsP3R activation, and vice versa, increasing Catot induced stable limit cycles even when there was no InsP3R activation in this system. As is evident from the lower panels in Figure 1, the frequency of the automaticity increased markedly with increasing Catot or pOInsP3R. The connection between the limit cycle identified in the bifurcation analysis and the cycle of spontaneous [Ca 2+ ] fluctuations in the bulk cytosol space (blk; [Ca 2+ ]blk) was examined by conducting numerical time integration of the same minimal model ( Figure 1C). Around the threshold Catot level  The connection between the limit cycle identified in the bifurcation analysis and the cycle of spontaneous [Ca 2+ ] fluctuations in the bulk cytosol space (blk; [Ca 2+ ] 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 2+ ] blk transient was firstly observed at 4.18 femtomole. The frequencies of the Ca 2+ transient increased with increasing Ca tot as shown in Figure 1A. These responses of the minimal model were completely reversible when Ca tot decreased.

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 Clcurrent (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].

Generation of APTDs and an Activation Threshold for Ca 2+ tot
The bifurcation analysis ( Figure 1) applied to the intracellular Ca 2+ dynamics in the absence of membrane ion fluxes revealed that spontaneous and repetitive Ca 2+ transients were initiated by increasing Catot above a certain threshold level (Catot,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 Catot was technically increased step by step by enlarging the scaling factor of ICabg. The spontaneous discharge of APTD appeared when Catot was increased above 3.878 femtomole in the absence of AR stimulation, whereas the threshold Catot decreased to 3.617 femtomole in the presence of 0.15 μM [ISO] and 0.15 μM [IP3]. Although the cycle length was markedly different, in both cases, the time course of Vm 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 2+ 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 2+ 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.

Generation of AP TD s and an Activation Threshold for Ca 2+ tot
The bifurcation analysis ( Figure 1) applied to the intracellular Ca 2+ dynamics in the absence of membrane ion fluxes revealed that spontaneous and repetitive Ca 2+ 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 3 ]. 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 2+ flux (J Ca ) from SR via RyR (passive leak component (J leak_SR ) and active release component (J rel_SR )), and via InsP 3 R (J InsP3R ) recovered in parallel with the replenishment of [Ca 2+ ] SRrl within the following 400-500 ms because of the Ca 2+ diffusion from SR uptake sites. The total continuous Ca 2+ efflux (J rel_SR ) finally evoked the final full Ca 2+ release as indicated by the rapid fall of [Ca 2+ ] SRrl through an increase in [Ca 2+ ] 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 2+ 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 2+ ] 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. The distribution of Catot 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 APTDs could be evoked by applying NA at 0.15 μM at the control sfICabg = 1.68. Figure 3C,D show APTD as in Figures 3A and 3B, and demonstrate that the interval between the two successive APTDs is much reduced (interval = ~500 ms), and the preceding TD evoked by the spontaneous Ca 2+ release is obvious before the rising phase of the AP. Again, no slow diastolic depolarization was observed, and a rapid decay in JCa via a cluster of RyRs (couplon) clearly precedes the onset of TD, and the second notch of JCa 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 sfIK1 by three-fold ( Figure 3D). The decrease in the TD amplitude in Figure 3D compared to that in Figure 3B   Alternatively, the AP TD s 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 TD s is much reduced (interval =~500 ms), and the preceding TD evoked by the spontaneous Ca 2+ 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 3 ] of 0.015 and 0.15 µM. The threshold Ca tot at the 0.15 µM [IP 3 ] was higher (3.86 and 3.83 femtomole) than the value of (3.65 and 3.59 femtomole) obtained at 0.015 µM [IP 3 ]. 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.

Separation of α1and β1-AR Influences on Ca tot under Resting Condition
The findings in both Figures 1 and 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 2+ 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 2+ from extracellular space to the SR space. Ca tot was markedly decreased by increasing the cytosolic [IP 3 ] ( Figure 4B), which was applied to represent the activation of α1-AR. The application of two concentrations of IP 3 caused a dose-dependent decrease in Ca tot as shown by removing the desensitization of the receptor at [IP 3 ] 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.

Separation of α1-and β1-AR Influences on Catot under Resting Condition
The findings in both Figures 1 and 3 indicate that the increase in Catot is the key factor for the initiation of the APTD 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 Catot at a saturating concentration of ISO (0.2 μM). In Figure 4A, the increase in Catot 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 Catot 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 Catot was soon converted to a decrease during the ISO application due to a gradual decrease in Catot through the NCX-mediated Ca 2+ 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 2+ from extracellular space to the SR space. Catot was markedly decreased by increasing the cytosolic [IP3] ( Figure 4B), which was applied to represent the activation of α1-AR.

Marked Latency before the Onset of Repetitive APTD 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 Catot induced by activating InsP3R (α1-AR) might counteract the β1-AR activation, which promotes the APTD burst by increasing Catot through ICaL and SERCA activation ( Figure 4A). This was the case, as shown in Figure 5A2, where the Catot rapidly decreased (by ~7%) after the start of AR stimulation. In this simulation, the control Catot level was set slightly lower than the threshold of initiating APTD (green vertical line) in the absence of AR stimulation. Factors involved in this response are shown in Figure  5A4. The open probability (pOInsP3R) of InsP3R (green) quickly increased from 0.00027 to a peak of

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 3 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 3 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 3 R can be partially activated by Ca 2+ . Figure 5(A3) indicates that the [Ca 2+ ] SRup and [Ca 2+ ] SRrl , the major store site of Ca tot , followed the similar time course as Ca tot . 0.00385 and then decayed to 0.00039 due to the desensitization of the α1-AR receptor activation. The Catot followed a time course similar to pOInsP3R ( Figure 5A2). The slightly higher desensitized pOInsP3R than the control was due to a rise in Catot. Note, the InsP3R can be partially activated by Ca 2+ . Figure  5A3 indicates that the [Ca 2+ ]SRup and [Ca 2+ ]SRrl, the major store site of Catot, followed the similar time course as Catot.
The time courses of the β1-AR stimulation of LCC (afCaL, red trace in Figure 5A4), SERCA (afSERCA, black), and Na + /K + pump (afNaK, blue) are biphasic because of delayed fractional (30%) desensitization of β1-AR. The discharge of APTD evoked a full CICR, and thereby [Ca 2+ ]SRrl decreased to a minimum level ( Figure 5B3). The [Ca 2+ ]SRup increased approximately parallel to the rise in [Ca 2+ ]blk, and showed a transient increase due to the uptake of Ca 2+ via SERCA at every Ca 2+ transient. The wash out of the influence of AR stimulant took an exponential time course and lasted for ~2 min.  Figure 5A4 are applicable to responses in (B,C). The rate of APTD discharge (/min/100) appears only in Figure 5B4. The horizontal green and brown lines are the threshold Catot measured in the absence and presence of AR stimulation, respectively.

Involvement of the Spontaneous Membrane Fluctuations in Determining the Latency
In the control run shown in Figure 5A, the PVC remained quiescent because the Catot level was still lower than the reference level (green line) of the spontaneous APTD discharge even after the desensitization of α1-AR. The AP burst was initiated only when the Catot was further increased by augmenting the β1-AR stimulation. In such case, a train of APTDs started at a given delay, different from the experimental finding of large variation in the 0-10 min latency [7]. For example, if scfICabg was increased from 1.68 to 2.2, Catot was still lower than the reference level, but the application of a saturating concentration of ISO (0.2 μM) evoked a train of APTD 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 Catot after the AR stimulation. In the The time courses of afI CaL (×0.5, red), af SERCA (black), af Na/K (blue), and pO InsP3R (×100, green) shown in Figure 5(A4) are applicable to responses in (B,C). The rate of AP TD discharge (/min/100) appears only in Figure 5(B4). The horizontal green and brown lines are the threshold Ca tot measured in the absence and presence of AR stimulation, respectively.
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 2+ ] SRrl decreased to a minimum level ( Figure 5(B3)). The [Ca 2+ ] SRup increased approximately parallel to the rise in [Ca 2+ ] blk , and showed a transient increase due to the uptake of Ca 2+ via SERCA at every Ca 2+ transient. The wash out of the influence of AR stimulant took an exponential time course and lasted for~2 min.

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 TD s 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 TD s. The size of the sporadic CICR was variable at individual CICR events, as evidenced by the miniature fluctuations in [Ca 2+ ] 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 TD s. 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 TD s 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 2+ 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 2+ 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 2+ 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 TD s were absent ( Figure 5A,C).

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 TD s 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 2+ influx due to the activation of I CaL and the uptake of Ca 2+ into SR by SERCA play a key role in determining the time-dependent increase in Ca tot . Note, the [Ca 2+ ] nd is largely dependent on [Ca 2+ ] SRrl when a continuous Ca 2+ leak is generated by both InsP 3 Rs and leak conductance of couplons.

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 2+ 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 Catot within the cell (Figure 1). The threshold level for the initiation of the repetitive APTD generation was slightly lower in the full model than in the model of cytosolic Ca 2+ 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 Catot, whereas the α1-AR activation temporary decreased Catot (Figure 4). The marked latency of several minutes for the start of repetitive APTD 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 2+ 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.

Co-Localization of InsP3R with RyRs in the Sub-Sarcolemmal Space Supporting Spontaneous CICR
The experimental studies in rat atrial myocytes demonstrated that a long lasting Ca 2+ release from SR is induced by IP3, and through the co-localization of InsP3R with RyRs, the transient Ca 2+ release is evoked [10]. Essentially the same mechanism was demonstrated when the InsP3R was activated by endothelin-1 (ET-1) [11]. The Ca 2+ extrusion through NCX is augmented by this transient increase in [Ca 2+ ] and thereby a temporal transient membrane depolarization is evoked due to the electrogenic stoichiomety of Na + /Ca 2+ exchange [19][20][21][22]. Mackenzie et al. directly recorded Ca 2+ sparks as well as premature APs in the presence of ET-1 at a high probability of 75% in trials in rat

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 2+ 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 2+ 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 2+ 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.
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 2+ as firstly described in the HuVEC model [16,47]. Ca 2+ accumulated in jnc gradually diffuses to iz and then to blk, in which myofilaments are located. Increased Ca 2+ 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 2+ uptake site of the network SR spread in the blk for computational simplicity. The major compartment of SR is SRup equipped with SERCA, and provides its extension (terminal cysterna) to form SRrl as the other counterpart of the dyadic junction. The SRrl site contains the Ca 2+ binding protein, calsequestrin, to increase the capacity of Ca 2+ release. The Ca 2+ released via InsP3R diffuses into the jnc. The ion channels and transporters are exposed to different [Ca 2+ ] after the Ca 2+ release, and their distribution to each compartment is given in the Appendix, but not shown in Figure 7. During the Ca 2+ release, [Ca 2+ ]SRrl is depleted to halt the Ca 2+ 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 INa, APTD is triggered. A series of AP discharge occurs when Catot within the cell is increased above a threshold level. In the present study, the extent of Ca 2+ overload was adjusted simply by magnifying the hypothetical ICabg assumed on the surface membrane. Simulating various mechanisms of the Ca 2+ 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 2+ concentration in the blk.

Co-Localization of InsP 3 R with RyRs in the Sub-Sarcolemmal Space Supporting Spontaneous CICR
The experimental studies in rat atrial myocytes demonstrated that a long lasting Ca 2+ release from SR is induced by IP 3 , and through the co-localization of InsP 3 R with RyRs, the transient Ca 2+ release is evoked [10]. Essentially the same mechanism was demonstrated when the InsP 3 R was activated by endothelin-1 (ET-1) [11]. The Ca 2+ extrusion through NCX is augmented by this transient increase in [Ca 2+ ] and thereby a temporal transient membrane depolarization is evoked due to the electrogenic stoichiomety of Na + /Ca 2+ exchange [19][20][21][22]. Mackenzie et al. directly recorded Ca 2+ 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 2+ sparks. The involvement of InsP 3 R was proven by recording larger amplitudes of premature AP when a membrane-permeable analogue of IP 3 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 2+ mechanisms suggested by the previous studies. In the present study, a functional coupling of InsP 3 R with RyRs was achieved by connecting the Ca 2+ flux through InsP 3 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 2+ ] within this jnc in the Hinch formalism of CaRU [14][15][16].

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 3 R at a much higher density than in the ventricular myocytes. InsP 3 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 3 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 3 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 3 R should be much higher in PVCs than in atrial myocytes. This well agrees with PVs showing high frequencies of AP TD s 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 3 R [10,11]. Note, the high frequency of AP TD generation is attributable to the frequency of the intracellular Ca 2+ 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 TD s. In conclusion, The AP TD s 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 TD s in rat.

Ca 2+ 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 2+ overload is also caused by the repetitive stimulation of ventricular myocytes in the presence of α1-NA stimulation, and the spontaneous Ca 2+ -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 2+ 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 α1plus β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 2+ release was evoked simply by increasing the Ca tot in the absence of the membrane ionic mechanisms.

The Latency and Frequency of AP TD s 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 2+ -mediated inactivation induced by a temporal release of SR Ca 2+ through InsP 3 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 2+ 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 2+ [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 3 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.

Coupling Several Layers of Physiological Mechanisms to Regulate Ca tot in PVCs
A regular time course of [Ca 2+ ] SRrl change (core cycle) is evoked when the cell is stimulated by an isolated electrical stimulus; the CICR transiently depletes [Ca 2+ ] SRrl at the terminal cysterna through excitation-contraction (EC) coupling. Then, the [Ca 2+ ] SRrl gradually recovers through the Ca 2+ diffusion from the Ca 2+ uptake site of SR. If Ca 2+ leaks through RyRs and InsP 3 R are added to this core cycle, any Ca 2+ accumulation in jnc potentially triggers the next event of Ca 2+ release via CaRU. The bifurcation analysis (Figure 1) indicated that a stable limit cycle of this Ca 2+ release is evoked by increasing the Ca tot beyond a certain level in the absence of Ca 2+ flux across the cell membrane.
If the membrane ionic fluxes are coupled with the core cycle of Ca 2+ oscillation, positive and negative feedback mechanisms are added. As a negative feedback, the transient Ca 2+ release via couplons primarily accelerates extrusion of Ca 2+ to the extracellular space via NCX. If an enlarged TD, evoked by the electrogenic Na + /Ca 2+ exchange, triggers an AP TD , the extra Ca 2+ influx through LCC activation can increase Ca tot , which may accelerate the core cycle of Ca 2+ 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 2+ dynamics. Basically, the simultaneous activation of both LCC and SERCA efficiently accumulates Ca tot through Ca 2+ uptake into SR. The activation of InsP 3 R through α1-AR activation strongly decreases Ca tot in combination with the Ca 2+ extrusion by NCX. This simple antagonistic relationship between α1and β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 2+ ] 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 α1and β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 3 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.

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 3 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]. 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 2+ 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 2+ domain called jnc, which allows a temporal accumulation of released Ca 2+ as firstly described in the HuVEC model [16,47]. Ca 2+ accumulated in jnc gradually diffuses to iz and then to blk, in which myofilaments are located. Increased Ca 2+ 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 2+ uptake site of the network SR spread in the blk for computational simplicity.

Intracellular Ca 2+ Compartments and Distribution of Ionic Channels and Transporters in the PVC Model
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 2+ binding protein, calsequestrin, to increase the capacity of Ca 2+ release. The Ca 2+ released via InsP 3 R diffuses into the jnc. The ion channels and transporters are exposed to different [Ca 2+ ] after the Ca 2+ release, and their distribution to each compartment is given in the Supplemental Materials, but not shown in Figure 7. During the Ca 2+ release, [Ca 2+ ] SRrl is depleted to halt the Ca 2+ 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 2+ overload was adjusted simply by magnifying the hypothetical I Cabg assumed on the surface membrane. Simulating various mechanisms of the Ca 2+ 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 2+ concentration in the blk.

Relationship Between Local [Ca 2+ ] nd and Spontaneous Ca 2+ Release
In the Hinch model, the [Ca 2+ ] nd is determined in different combinations of LCC and couplons under the assumption of instantaneous Ca 2+ distribution [14,15].
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.
where G L and G R are Ca 2+ conductivity through the LCC and couplons, respectively; and G diff represents the conductivity of Ca 2+ 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 2+ ] 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).
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 2+ fluxes in the jnc space, Ca 2+ release through couplon and InsP 3 R, and Ca 2+ diffusion between the jnc and iz, respectively. Under a Ca 2+ -overload condition, the J Ca_rel through the basal openings of the couplon or activated InsP 3 R increases in proportion to the level of [Ca 2+ ] SRrl above the threshold level. The amplitude of Ca 2+ flux is determined by Equations (5) and (6).
For Ca-buffers in jnc, iz, and blk, see Supplemental Materials.

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 2+ compartment (Table 1). 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 2+ -mediated inactivation of LCC is nd and the inactivation rate (k oc ) is described as: 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 2+ ]-dependence in nd or other [Ca 2+ ]. 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 2+ 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 1 , C 2 , 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).

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 α1and β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 3 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.

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): The time course of Fr PKA after the onset of β1-AR stimulation was calculated using Equation (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).
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. a f (t) = base + delta·Fr PKA(t) .

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.
To satisfy the thermodynamic constraint of: 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 I NaK ,

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 2+ 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 2+ flux into SR was calculated by Equation (17). The same time constant τ = 1 s was used for the β1-AR regulation of Na + /K + .
The whole cell J SERCA was determined as a sum of control J SERCA_c and activated J SERCA_a .

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)).

α1-Adrenergic Signal, [IP 3 ]
We failed to find an α1-AR model to incorporate in the PVC model. In this study, the timedependent activation and desensitization of the receptor was simply represented by changes in the active second messenger [IP 3a ] using a concentration [IP 3 ] (t) , and an inactivation (i) parameter of the α1-AR at time t. To obtain [IP 3 ], equal to the control 0.015 mM after the desensitization, a lower limit of 0.1 was applied to the range of i ∞ : The [IP 3 ] (t) at time t is also given by an exponential function: 4.4.6. InsP 3 R In the cardiac myocytes, the type II homologous isomer is the major component of InsP 3 R [59] and is expressed on the SR membrane near the T-tubules [7]. Thus, in the model, the InsP 3 R was exposed to jnc, so that its Ca 2+ release plays a critical role in determining the bias level for the CICR. The open probability pO InsP3R was calculated using the InsP 3 R-II model [60]. The influence of Ca 2+ on the pO InsP3R was much less than that of IP 3 over the range of [Ca 2+ ] during the Ca 2+ transient. The resting level of [IP 3 ] was set to 0.015 µM according to Cooling et al. [17] so that J InsP3R was virtually zero. The [IP 3 ] was increased 10-fold (0.15 µM) during the stimulation of α1-AR. The Ca 2+ flux via InsP 3 R channel was calculated using Equation (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.

Simulation of the Random Events of CICR
If [Ca 2+ ] 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 2+ ]/dt, the average time course of [Ca 2+ ] 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 2+ ] 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 2+ ] jnc is only used to determine the opening rate of a couplon, but does not interfere with the mass conservation of Ca 2+ dynamics.