Stretch-activated current in human atrial myocytes and Na+ current and mechano-gated channels’ current in myofibroblasts alter myocyte mechanical behavior: a computational study

Background The activation of stretch-activated channels (SACs) in cardiac myocytes, which changes the phases of action potential repolarization, is proven to be highly efficient for the conversion of atrial fibrillation. The expression of Na+ current in myofibroblasts (Mfbs) regenerates myocytes’ action potentials, suggesting that Mfbs play an active role in triggering cardiac rhythm disturbances. Moreover, the excitation of mechano-gated channels (MGCs) in Mfbs depolarizes their membrane potential and contributes to the increased risk of post-infarct arrhythmia. Although these electrophysiological mechanisms have been largely known, the roles of these currents in cardiac mechanics are still debated. In this study, we aimed to investigate the mechanical influence of these currents via mathematical modeling. A novel mathematical model was developed by integrating models of human atrial myocyte (including the stretch-activated current, Ca2+–force relation, and mechanical behavior of a single segment) and Mfb (including our formulation of Na+ current and mechano-gated channels’ current). The effects of the changes in basic cycle length, number of coupled Mfbs and intercellular coupling conductance on myocyte mechanical properties were compared. Results Our results indicated that these three currents significantly regulated myocyte mechanical parameters. In isosarcometric contraction, these currents increased segment force by 13.8–36.6% and dropped element length by 12.1–31.5%. In isotonic contraction, there are 2.7–5.9% growth and 0.9–24% reduction. Effects of these currents on the extremum of myocyte mechanical parameters become more significant with the increase of basic cycle length, number of coupled Mfbs and intercellular coupling conductance. Conclusions The results demonstrated that stretch-activated current in myocytes and Na+ current and mechano-gated channels’ current in Mfbs significantly influenced myocyte mechanical behavior and should be considered in future cardiac mechanical mathematical modeling.

Our previous study has found that I Na_myofb and I MGC_Mfb regenerated APs in myocytes and Mfbs [28]. In this study, we aimed to investigate the role of I SAC , I Na_Mfb , and I MGC_Mfb in the mechanical contraction of cardiac myocyte. Simulation results of human atrial myocyte segment mechanical dynamics with different gap-junctional conductance (G gap ), number of coupled Mfbs and basic cycle lengths (BCLs) were examined.

Results
Effects of I SAC , I Na_Mfb , and I MGC_Mfb on atrial myocyte AP, [Ca 2+ ] i , and the normalized force Figure 1 shows the combinational effects in five groups (see "Simulation protocol" section) of I SAC , I Na_Mfb , and I MGC_Mfb on the membrane potential, intracellular Ca 2+ concentration and the normalized force (F norm ) of myocytes with a G gap of 3 nS and a BCL of 1 s. For myocytes, coupling Mfbs (Group 2-5) resulted in gradual decrease of myocyte membrane potential amplitude (V max ) and APD at 90% repolarization (APD 90 ), and increase of the resting myocyte membrane potential (V rest ) depolarization (Fig. 1a). Meanwhile, a spontaneous excitement was emerged in Group 5, in which the peak [Ca 2+ ] i dropped significantly (Fig. 1b), indicating that I SAC , I Na_Mfb , and I MGC_Mfb could result in discordant alternans. From the traces of F norm (Fig. 1c), it can be observed that the peak F norm increased after myocyte coupled to Mfbs. It was increased by 7.6% (Group 2), 14.5% (Group 3), 38.7% (Group 4), and 19.2% (Group 5) as compared to the control (Group 1). It was remarkable that Group 4 got the biggest F norm increment, which meant F norm of myocytes could be significantly enhanced by the combination of I SAC and I Na_Mfb . However, the increment relatively declined in 10 10  Effects of I SAC , I Na_Mfb , and I MGC_Mfb on atrial myocyte segment mechanical parameters Traces of F SE , F PE , F segment , l CE , l PE , and l SE obtained in five groups for the simulations of isosarcometric contraction with sarcomere length of 1.78 µm are displayed in Fig. 2, and the ones for simulations of isotonic contraction with applied force of 10 mN per square millimeter (mN/mm 2 ) are displayed in Fig. 3. In isosarcometric contraction (Fig. 2), F PE and l PE in five groups were constant. Peak F SE increased when myocyte coupled to Mfbs. It was increased by 7.2% (Group 2), 13.8% (Group 3), 36.6% (Group 4), and 18.5% (Group 5) as compared to the control (Group 1). Using Eq. (14) (see "Mechanical behavior of a single segment" section), the increments in F segment were the same as those in F SE . On the contrary, l CE decreased when myocyte coupled to Mfbs. The minimum of l CE was dropped by 6.4% (Group 2), 12.1% (Group 3), 31.5% (Group 4), and 16.2% (Group 5) as compared to the control (Group 1). Using Eq. (15), the changes in l CE were equal to those in l SE . Similar to Fig. 1, Group 4 had the most significant change, indicating that the combination of I SAC and I Na_Mfb played a significant role in determining myocyte segment mechanical behavior.
In isotonic contraction (Fig. 3), F segment in five groups was constant. Like Fig. 2a, the peak value of F SE also increased when Mfbs were coupled. It was increased by 1.6% (Group 2), 2.7% (Group 3), 5.9% (Group 4), and 2.5% (Group 5) as compared to the control (Group 1). The increments were smaller than those in isosarcometric contraction. According to Eq. (14), the decrements in F PE were equal to the increments in F SE . The minimum of l CE was dropped by 4.3% (Group 2), 7.4% (Group 3), 15.8% (Group 4), and 6.3% (Group 5) as compared to the control (Group 1), while the peak value of l SE was increased by 0.6%, 1.4%, 3.7%, and 0.9%, and finally, l PE were declined by 6.8%, 11.4%, 24.0%, and 10.0%. I SAC , I Na_Mfb , and I MGC_Mfb not only changed the extreme values of l CE , l PE , and l SE , but also altered l CE and l SE in the resting stage. For example, the value of l SE in resting period was 0.0021 μm (Group 1), 0.0013 μm (Group 2), 0.0021 μm (Group 3), 0.0051 μm (Group 4), and 0.0032 μm (Group 5), respectively. Similarly, Group 4 has the most impact on the mechanical parameters.
To investigate the effects of I SAC , I Na_Mfb , and I MGC_Mfb on the extreme values of atrial myocyte segment mechanical parameters, we simulated five groups with different BCLs, Mfb-M ratios, and G gap in both isosarcometric contraction (Fig. 4) and isotonic contraction (Fig. 5).
From the traces of the extremum of F SE , F segment , l CE , and l SE for BCL = 0. , and l SE [MAX] in Group 4 reached their maximums, and l SE[MAX] obtained its minimum. This suggests that I SAC together with I Na_Mfb had the key influence on myocyte mechanical parameters. The influence disappeared in Group 5, suggesting that the role of I MGC_Mfb in myocytes was opposite compared to I SAC . As BCL longer than 1 s, each parameter in Group 1 to 3 has increased or decreased as a same trend, whereas fluctuated in Group 4 and 5. These phenomena might be attributed to that I SAC and I Na_Mfb enhanced atrial myocytes excitability and triggered spontaneous excitements at large BCLs. [Ca 2+ ] i , the vehicle of EMC, also fluctuated, driving the undulation of mechanical behavior. Our results demonstrated that introducing currents through SACs in myocytes and currents through MGCs in Mfbs in cardiac modeling could lead to different simulation results. In fact, the stretch ability and contractility of myocytes in fibrotic heart were quite different from those in normal heart. Integrated I SAC and I MGC_Mfb in cardiac simulation could help obtain more accurate and closer to experimental results. Figure 4i-l shows the extremum of four parameters with G gap ranging from 0.5 to 8 nS. Parameters in Group 1 were also constants. The variance of five groups was less than those in other settings ( Fig. 4a-h), suggesting the relative small effects of G gap on myocyte mechanical parameters. The traces of Group 2 to 5 were mostly distributed over one side of Group 1. Parameters at each G gap in Group 4 got the highest or lowest value among five groups, and parameters in Group 5 took the second place, indicating that I SAC , I Na_Mfb , and I MGC_Mfb played a strong role in atrial myocyte mechanical behavior. The extremum of F SE , F PE , l CE , l SE , and l PE as functions of BCL, Mfb-M ratio, and G gap in isotonic contraction are showed in Fig. 5.
In Fig. 5a-e, the parameters among five groups were close to each other. As BCL increased, the values in Group 1, 4 and 5 first increased and then decreased. In these groups, pure myocyte or integrating I SAC and I MGC_Mfb in fibrotic myocyte were more likely to cause discordant alternans and mechanical parameters fluctuation at big BCLs.
In Fig. 5f-j, the parameters in Group 2 always stayed over one side of Group 1 as the coupled ratio increased, while the parameters in other groups finally converged over the other side. It suggested that integrating I SAC , I Na_Mfb , and I MGC_Mfb in fibrotic myocyte significantly influenced the myocyte segment mechanical behavior at large coupled ratios. In Fig. 5k-o, the parameters in Group 1 were constant, and Group 2, 3, and 5 had the similar traces, while Group 4 got the highest or lowest values. Therefore, I SAC , I Na_Mfb , and I MGC_Mfb had relative influences on mechanical parameters at large G gap .

Discussion
This study investigated the roles of I SAC , I Na_Mfb , and I MGC_Mfb in myocyte segment mechanical behavior. To address these issues, computational simulations of the coupled Mfb-M system were performed by employing a combination of models of the human atrial myocyte (including I SAC ) and Mfb (including I Na_Mfb and I MGC_Mfb ), as well as models of Ca 2+ -force relation and myocyte mechanical segment. Specifically, effects of these currents with changes in (1) BCL, (2) the number of coupled Mfbs, and (3) G gap on atrial myocyte segment mechanical parameters were investigated. The integration of I SAC , I Na_ Mfb , and I MGC_Mfb could result in (1) decreased V max and APD 90 , increased V rest depolarization, and spontaneous excitements even discordant alternans at large BCLs, and (2) increased peak value of F SE , F segment , and l SE and decreased valley value of l CE in isosarcometric contraction, and increased peak value of F SE and l SE and decreased valley value of F PE , l CE , and l PE in isotonic contraction. Moreover, I SAC and I MGC_Mfb have relative effects on myocyte segment mechanical parameters.

Effects of I SAC , I Na_Mfb , and I MGC_Mfb on atrial myocyte segment mechanical properties
Effects of I SAC , I Na_Mfb , and I MGC_Mfb on the excitability of human atrial myocytes have been discussed in our previous study [28]. Here, we discussed the roles of these currents in myocyte segment mechanical behavior.
EMC and MEF were two known effects [7], but the physical role of MEF in EMC was still poorly understood. In general, I SAC , handling as the major mechanisms of the MEF, was reported to enhance the early phase of AP repolarization and prolong or delay the final phase of repolarization [9,35,36]. But the impact of I SAC on cardiac mechanics, to our best knowledge, has been rarely studied so far. In our present study, the stretchactivated currents had the most significant influence on myocyte segment mechanical parameters in both isosarcometric contraction and isotonic contraction.
For cardiac Mfbs, many studies have verified that mechanical cues activated cardiac Mfbs and led to increased production of extracellular matrix [37,38]. Mfbs were regarded as a critical determinant of cardiac mechanics. Previous studies have used computational modeling to demonstrate the acute mechanical effects on cardiac fibroblast structure and organization [39,40]. They found that an axial strain environment could guide fibroblast proliferation, orientation, and migration [31,41,42]. Several groups have simulated cell compaction of collagen gels by calculating mechanical equilibrium between each cell's contractile forces and nearby collagen fibers' mechanical properties. They reported that cellular organization is tightly linked to the mechanical feedback loop between cells and matrix [29,30]. These studies were all about the stretch-induced responses of quiescent cardiac Mfbs. However, the inverse process, i.e., the Mfbs-induced responses of cardiac mechanics, has not been widely explored. Our results showed that coupling Mfbs changed myocytes mechanical properties. In addition, we compared the results of before and after adding I Na_Mfb and I MGC_Mfb in the Mfb model, and found that the effects of I MGC_Mfb on the force of atrial myocytes were contrary to I SAC . For I Na_Mfb , many studies have been conducted to investigate how this current could influence Mfb proliferation [18,43]. Our results showed that I Na_Mfb decreased V max and APD 90 and increased V rest depolarization in myocytes. This depolarization changed diastolic Ca 2+ levels and then altered myocytes mechanical behavior.
For I MGC_Mfb , experimental data have indicated that cardiac fibroblasts expressed functional MCGs, contributing to the cardiac MEF both under physiological and pathophysiological conditions [44,45]. We assumed that it could affect myocytes mechanical characteristics like I SAC . Our results supported this hypothesis. In our simulations, I MGC_Mfb altered myocytes mechanical behavior. Interestingly, the effects of I MGC_Mfb and I SAC on myocyte segment mechanical parameters seemed to be opposite. Myocytes stretch activated I SAC and enhanced the influence on mechanical parameters, while Mfbs compression activated I MGC_Mfb and weakened the influence. Moreover, MGCs were activated by fibroblast compression and inactivated by fibroblast stretch [21], implying that I MGC_Mfb should be integrated in cell modeling only during cell compression, such as fibroblasts/Mfbs compression caused by stretching and dilatation of surrounding cardiac myocytes.
Mfb was a critical determinant of cardiac mechanics. Previous studies have demonstrated that abnormal quantity or organization of Mfb could lead to both systolic and diastolic dysfunction [12,30,46]. Besides, previous modeling work suggested that Mfb-M coupling contributed to arrhythmia formation [25,47]. The key factors included BCL, the number of coupled Mfbs, and G gap . Here, we integrated I Na_Mfb , I SAC , and I MGC_ Mfb into Mfb-M coupling and compared their effects on myocyte mechanical properties in different settings of BCL, Mfb-M ratio, and G gap . To the best of our knowledge, this has not been examined before. With BCL, Mfb-M ratio, and G gap increasing, impacts of these currents on the extremum of myocyte mechanical parameters became greater, as summarized in Figs. 4 and 5.

Limitations
Two limitations in the present study should be mentioned. First, functional roles of SACs in Mfbs were not considered. Direct proof of mechanoactivation of mechanosensitive channels in cardiac Mfbs was limited. A handful of experimental studies have found that mechanical cues could lead to the opening of so-called SACs, and transient receptor potential canonical channels were candidates for the stretch-activated currents measured in cardiac fibroblasts [48,49]. However, the current-voltage relation of I SAC in Mfbs needs further study. Second, the breadth of this computational study needs to be extended. Our work focused on the scale of local cell-cell interactions. Other scales, such as scales of subcellular signaling, cell-matrix interactions, tissue remodeling, and organ level conduction properties, were not included in this preliminary study. In fact, processes across these scales did not occur in isolation but operated as an interconnected system with every level passing information to other levels. Therefore, multi-scale modeling frameworks still need to be developed, although they brought computational challenges, and such models involving cardiac Mfbs and fibrosis were still rare.

Conclusions
This study demonstrated the combinational effects of I SAC in myocytes and I Na_Mfb and I MGC_Mfb in Mfbs on myocyte mechanical properties. Our results showed that the addition of I SAC , I Na_Mfb , and I MGC_Mfb regulated the peak and valley values of myocyte mechanical parameters in both isosarcometric contraction and isotonic contraction. Effects of these currents on the extremum of myocyte mechanical parameters become more evident as BCL, Mfb-M ratio, and G gap increased. The effects proved that the stretch-activated current in atrial myocyte and Na + current and mechano-gated channels in Mfbs should be considered in future pathological cardiac mechanical mathematical modeling, such as atrial fibrillation and cardiac fibrosis.
In following sections, the details of each component of the model will be described.

The model of Mfb-M coupling
The Mfb-M coupling will be modeled based on [26], with the differential equations for the membrane potential of cardiac Mfb and myocyte are given by where V Mfb,i and V M represent the transmembrane potential of the ith coupled Mfb and the human atrial myocyte, C m,Mfb and C m,M represent the membrane capacitance of the Mfb and the myocyte, I Mfb,i and I M represent the transmembrane current of the ith coupled Mfb and the human atrial myocyte, and G gap represents the gap-junctional conductance. It is also noted that a negative I gap [i.e., G gap (V Mfb,i -V M )] indicates that the current is flowing from the myocyte into the ith Mfb, and n is the total number of coupled Mfbs.

Mathematical model of the human atrial myocyte
The mathematical model of the human atrial myocyte developed by Maleckar et al. [50], which is based on experimental data and has correctly replicated APD restitution of the adult human atrial myocyte, was adopted in this study. To examine the influence of the stretch on myocyte AP, the original model from Maleckar et al. is modified with the total ionic current of myocyte (I M ) given as where I Na is fast inward Na + current, I CaL L-type Ca 2+ current, I t transient outward K + current, I Kur sustained outward K + current, I K1 inward-rectifying K + current, I K,r rapid delayed rectifier K + current, I K,s slow delayed rectifier K + current, I B,Na background Na + current, I B,Ca background Ca 2+ current, I NaK Na + -K + pump current, I CaP sarcolemmal Ca 2+ pump current, I NaCa Na + -Ca 2+ exchange current, I SAC stretch-activated current, and I Stim stimulated current.

The model of I SAC
Kuijpers et al. [55] have conducted experimental studies and reported that I SAC in atrial myocytes is permeable to Na + , K + , and Ca 2+ [33], and defined as where I SAC,Na , I SAC,K , and I SAC,Ca represent the contributions of Na + , K + , and Ca 2+ to I SAC , respectively. These currents are defined by the constant-field Goldman-Hodgkin-Katz current equation [56]. where F is Faraday's constant, Vol i cytosolic volume, I di Ca 2+ diffusion current from the diffusion-restricted subsarcolemmal space to the cytosol, I up sarcoplasmic reticulum Ca 2+ uptake current, I rel sarcoplasmic reticulum Ca 2+ release current, and O buffer occupancy.

The model of the Ca 2+ -force relation
The model 4 of isometric force generation in cardiac myofilaments proposed by Rice et al. was adopted to model the Ca 2+ -force relation [51,52] represents the total troponin low-affinity site concentration, and k + ltrpn and k − ltrpn are the Ca 2+ on-and off-rates for troponin low-affinity sites.

The model of human atrial Mfb
The electrophysiological model of human atrial Mfb proposed by MacCannell et al. [26] was used in the present study. It includes time-and voltage-dependent K + current (I Kv_Mfb ), inward-rectifying K + current (I K1_Mfb ), Na + -K + pump current (I NaK_Mfb ), and Na + background current (I B,Na_Mfb ).
In addition, I Na_Mfb and I MGC_Mfb are added in the Mfb model. According to our previous work [28], equations of I Na_Mfb and I MGC_Mfb are formulated as where g Na,Mfb is the maximum conductance of I Na_Mfb (0.756 nS), E Na,Mfb the Nernst potential for Na + ions, [Na + ] c,Mfb the Mfb extracellular Na + concentration (130.011 mM), [Na + ] i,Mfb the Mfb intracellular Na + concentration (the initial value is set as 8.5547 mM), and m Mfb and j Mfb the activation and inactivation parameters, respectively. To follow the experiment data [18], j has been modified as j 0.12 . g MGC,Mfb is the maximum conductance of I MGC_Mfb (0.043 nS), and E MGC,Mfb is the reversal potential of MGCs (selected a value close to 0 mV) [21].

Mechanical behavior of a single segment
The mechanical behavior of a single segment in our model is based on the classical threeelement rheological scheme [53,54].
As shown in Fig. 6, active force (F CE ) is generated by the contractile element (CE), and passive forces (F SE , F PE ) are generated in a serial elastic element (SE) and a parallel elastic element (PE). F segment is the total force generated by the segment. The element lengths are l CE , l SE , and l PE . During mechanical equilibrium, F CE , F segment , and l PE are defined as (10)

Simulation protocol
We performed single-cell simulations with constant sarcomere length (isosarcometric contraction) and constant applied force (isotonic contraction) to investigate the effects of I SAC , I Na_Mfb , and I MGC_Mfb on myocyte mechanical properties. Five groups were simulated sequentially: one atrial myocyte without Mfb coupling (Group 1), one atrial myocyte coupled to two Mfbs without I SAC , I Na_Mfb , and I MGC_Mfb (Group 2), one atrial myocyte coupled to two Mfbs with I Na_Mfb (Group 3), one atrial myocyte coupled to two Mfbs with I SAC and I Na_Mfb (Group 4), and one atrial myocyte coupled to two Mfbs with I SAC , I Na_Mfb , and I MGC_Mfb (Group 5).
First, simulations were carried out at a constant G gap of 3 nS and a BCL of 1 s. Thereafter, the coupled system was paced with (1) BCLs from 0.1 to 2 s, (2) G gap from 0.5 to 8 nS, and (3) number of coupled Mfbs from 1 to 8, to investigate the role of BCL, G gap , and Mfbs in myocyte mechanical parameters. The maximum or minimum of F SE , F PE , F segment , l CE , l SE , and l PE at different BCL, G gap , and Mfbs number were examined.
To ensure the coupled system reached steady-state, stimulation was repeated for 20 cycles. Results from the last cycle in each simulation were used for subsequent analyses. All state variables of the coupled model were updated by means of the forward Euler method. The time step was set to be 10 μs to ensure numerical accuracy and stability. More information on "Methods" is available in the Additional file 1.