Similar articles

Abstract: In this study, the effects of cardiac fibroblast proliferation on cardiac electric excitation conduction and mechanical contraction were investigated using a proposed integrated myocardial-fibroblastic electromechanical model. At the cellular level, models of the human ventricular myocyte and fibroblast were modified to incorporate a model of cardiac mechanical contraction and cooperativity mechanisms. Cellular electromechanical coupling was realized with a calcium buffer. At the tissue level, electrical excitation conduction was coupled to an elastic mechanics model in which the finite difference method (FDM) was used to solve electrical excitation equations, and the finite element method (FEM) was used to solve mechanics equations. The electromechanical properties of the proposed integrated model were investigated in one or two dimensions under normal and ischemic pathological conditions. fibroblast proliferation slowed wave propagation, induced a conduction block, decreased strains in the fibroblast proliferous tissue, and increased dispersions in depolarization, repolarization, and action potential duration (APD). It also distorted the wave-front, leading to the initiation and maintenance of re-entry, and resulted in a sustained contraction in the proliferous areas. This study demonstrated the important role that fibroblast proliferation plays in modulating cardiac electromechanical behaviour and which should be considered in planning future heart-modeling studies.

成纤维细胞增殖改变心脏兴奋传导和收缩的仿真研究

研究目的：通过建立心肌-成纤维细胞电力耦合模型，探讨成纤维细胞增殖对心脏兴奋传导和力学收缩的影响。
创新要点：基于人心室肌细胞的实验数据，修正了ten Tusscher等人发表的心室肌模型，构建心肌细胞-成纤维细胞的耦合模型，仿真了成纤维细胞对心肌组织电生理及力学的影响。通过改变耦合模型中重要参数的数值，如网格分辨率、成纤维细胞模型参数、缝隙连接电导等，观察其对心脏去极化、复极化和动作电位周期的影响。
重要结论：在细胞水平上，耦合成纤维细胞使心室肌细胞的动作电位周期延长，主动张力峰值下降（见图4、5）。在组织水平上，成纤维细胞增殖降低兴奋波传导速度并引发传导阻滞（见图6），降低成纤维增殖区域的应变，延长组织的去极化和复极化（见图8），并维持折返（见图10）。

关键词：心脏模型；电力耦合；成纤维细胞增殖

1.  Introduction

 To explore cardiac dynamics in both normal and pathological situations, cardiac modeling, including electrophysiology and mechanics, has been used as an important alternative to experimental studies. Many cardiac models, from cellular level models based on ion channels to tissue level models for ventricles, atria and the whole heart, have been proposed to study the electrophysiological and mechanical activities of the heart (Cherry and Fenton, 2011; Clayton et al., 2011; Trayanova, 2011; Trayanova et al., 2011; Carusi et al., 2012).

 Owing to the close relationship between myocardium electrical and mechanical activities, cardiac-coupled electromechanical models were developed, in from one to three dimensions, for both normal and pathological situations (Usyk and McCulloch, 2003a; 2003b; Nash and Panfilov, 2004; Nickerson et al., 2005; Niederer and Smith, 2008; 2009; Kerckhoffs et al., 2010; Niederer et al., 2011; Trayanova, 2011; Trayanova et al., 2011; Kuijpers et al., 2012). These coupled models fall mainly into two categories, namely, weak and strong coupling models (Kerckhoffs et al., 2006). The weak coupling method assumes that the coupling effects between deformation and electrophysiology are nominal, while the strong coupling method calculates the action potential (AP), active tension (T _a), and deformation simultaneously (Niederer and Smith, 2008). In terms of the weak coupling method proposed by Kerckhoffs et al. (2006), Xia et al. (2005) proposed an electromechanical biventricular model to analyze ventricular wall motion. The corresponding changes in the electrocardiogram (ECG) were simulated by Xia et al. (2006) and in the magnetocardiogram (MCG) by Shou et al. (2011). In general, the strong coupling assumption is more consistent with experimental studies, but has decreased numerical stability and increased computational load.

 However, these cardiac models describe only the properties of cardiac myocytes. Other types of cells, such as those in vasculature and connective tissues, (Manabe et al., 2002), are not considered. Recent studies have revealed that fibroblasts, as one kind of connective tissue, are the most abundant cell type by number in cardiac tissues, comprising about 70% of all cells (Banerjee et al., 2007). Because of their small size, they make up only about 6% of the total cardiac volume (Rossi, 2001). They affect the restitution properties of cardiac tissues and spiral wave stability (de Bakker and van Rijen, 2006; Spach et al., 2007; ten Tusscher and Panfilov, 2007; Xie et al., 2009a), especially during the process of ageing and in various cardiac diseases. Once the percentage of fibroblasts in the heart increases up to 10%–35%, a remodeling of the cardiac structure occurs, increasing muscle stiffness and reducing the coupling between adjacent muscle fiber bundles (Rossi, 2001; ten Tusscher and Panfilov, 2007). This process is strongly correlated with atrial and ventricular tachyarrhythmias and sudden cardiac death (John et al., 2004; Assomull et al., 2006; Pellman et al., 2010). The elastic modulus rises from (18±2) kPa in the normal heart muscle of non-infarcted animals to (55±15) kPa in significant fibroblast proliferous muscle of infarcted animals, an almost three-fold increase (Berry et al., 2006). Therefore, the functional roles of the fibroblasts in cardiac electrical and mechanical activities have attracted increasing interest, and have been explored in experimental studies (Miragoli et al., 2007; Zlochiver et al., 2008) and computational models (MacCannell et al., 2007; Tanaka et al., 2007; Jacquemet and Henriquez, 2008; Zlochiver et al., 2008; Sachse et al., 2009). Both experimental and computational studies (Miragoli et al., 2006; Jacquemet and Henriquez, 2008; Zlochiver et al., 2008) have illustrated that non-monotonic changes in conduction velocity (CV) are caused by increasing the fibroblast population. ten Tusscher and Panfilov (2007) showed that CV decreased from 0.7 to 0.3 m/s as the percentage of fibroblast proliferous cells rose from 0 to 40% in two-dimensional (2D) tissue. In simulation studies, it was reported that myocyte-fibroblast coupling modulated AP morphology (MacCannell et al., 2007; Jacquemet and Henriquez, 2008; Sachse et al., 2009). Coupling two or four fibroblasts to a ventricular myocyte resulted in a decrease in action potential duration (APD) from a control value of 263 ms to 195 or 155 ms, respectively (MacCannell et al., 2007). These studies suggest that the effects of fibroblasts on cardiac electrophysiology and mechanics should not be ignored. So far, two electrophysiological models of ventricular fibroblasts, a passive model (Xie et al., 2009b) and an active model (MacCannell et al., 2007), have been developed. However, to the best of our knowledge, there is no mechanical model for cardiac fibroblasts, although some experimental studies on the mechanics of the fibroblasts have been reported (Brown et al., 1998; Eastwood et al., 1998). These studies indicated that fibroblasts have the ability to maintain a tensional homeostasis of about (40–60)×10⁻⁵ N per million cells.

 In this study, a strongly coupled myocardial-fibroblastic electromechanical model is proposed to investigate the effect of fibroblasts on cardiac excitation conduction and contraction. The proposed model integrates the coupling between a fibroblast model (Xie et al., 2009a), a myocardial electrophysiological model (ten Tusscher et al., 2004), and a mechanical model (Rice et al., 1999) at the cellular level, with excitation conduction and elastic mechanics (Nash and Panfilov, 2004) at the tissue level. To simulate cellular electromechanical coupling, equations for calcium buffering developed by Iribe et al. (2006) were incorporated into the ten Tusscher model (ten Tusscher et al., 2004). The developed model was first validated by comparing the simulated results of the excitation conduction and mechanical properties to those from previous electromechanical coupling models. Using our model, we investigated the influence of diffuse fibroblast proliferation on wave propagation. Then, we analyzed the importance of various parameters, including tissue sizes, fibroblast areas, mesh resolution, parameters in the cellular fibroblast model, and the gap junctional coupling between the myocytes/fibroblasts, in determining the spatial distribution of depolarization, repolarization and APD. Strain maps and temporal traces of strain at different points of both normal and fibroblast proliferous tissues are illustrated. Electromechanical models of a central point stimulus and re-entry were simulated to investigate the effects of the fibroblasts on cardiac excitation conduction and mechanical contraction in normal and ischemic conditions.

2.  Materials and methods

2.1.  Model framework

 The detailed framework of the proposed strongly coupled myocardial-fibroblastic electromechanical model comprises two parts: the electrophysiological and the mechanical models of myocytes and fibroblasts at the cellular level, in addition to models of the excitation conduction and contraction of the fibroblast-myocyte coupling at the tissue level (Fig. 1). The cell model stems from the ten Tusscher electrophysiological model (ten Tusscher et al., 2004), the Rice model (Rice et al., 1999), and a passive fibroblast model (Xie et al., 2009a). The tissue model includes cardiac excitation conduction and finite deformation, as described by Nash and Panfilov (2004). Details of the individual components and the modification of each model are described as follows.

2.1.1.  Myocyte-fibroblast electrophysiological modeling at the cellular level

 For myocyte electrophysiology, the ten Tusscher model was employed in this study, as it accurately represents the detailed characteristics of the human ventricular cell (ten Tusscher et al., 2004). The only modification made was to the Ca²⁺ handling. To couple the electrophysiology with the mechanics of the cardiac myocyte, additional Ca²⁺ buffers were integrated into the Ca²⁺ dynamics to achieve an electromechanical coupling at the cellular level. The equations for the Ca²⁺ buffers are expressed as (Iribe et al., 2006) (1) , (2) , where Cmdn_Ca and Cmdn_tot are the individual and total concentrations, respectively, of Ca²⁺ bound to calmodulin; α _cmdn and β _cmdn are the Ca²⁺ on-rate and off-rate, respectively, for calmodulin; Ca_i is the cytosolic Ca²⁺ concentration; Trpn_Ca and Trpn_tot are the individual and total concentrations, respectively, of Ca²⁺ bound to troponin; α _trpn and β _trpn are the Ca²⁺ on-rate and off-rate, respectively, for troponin; and Force_norm is the normalized force.

 Correspondingly, the cytosolic concentration of Ca²⁺ in the ten Tusscher model was modified as (3) , where I _leak is a leakage flux from the sarcoplasmic reticulum to the cytoplasm; I _up is the sarcoplasmic reticulum (SR) Ca pump taking up calcium in the sarcoplasmic reticulum; I _rel is the calcium-induced calcium release (CICR) flux; I _CaL is the L-type calcium flux; I _bCa is the background calcium flux; I _pCa is the Ca²⁺ pump flux; I _NaCa is the Na⁺/Ca²⁺ exchanger flux; V _C is the cytoplasmic volume; and, F is the Faraday constant. More details about other parameters and equations in the ten Tusscher model can be found in ten Tusscher et al. (2004) and Iribe et al. (2006).

 For fibroblast electrophysiology, two previously proposed models can simulate the main characteristics of the fibroblasts (Chilton et al., 2005; Shibukawa et al., 2005; Xie et al., 2009b). The main difference between the models is in their representation of I _fb (fibroblast ionic current). The passive model assumes I _fb is one current (Xie et al., 2009a), while the active model describes it as a composite of several ionic currents (Chilton et al., 2005; Shibukawa et al., 2005). In this study, the passive model was applied (Xie et al., 2009a): (4) , where G _f is the membrane conductance, V _f is the membrane potential of the fibroblast, and E _f is the reversal potential. The G _f ranges from 0.1 to 4 nS (Kohl et al., 1994) and E _f from −60 to 0 mV (Kamkin et al., 1999; Kohl, 2003; Camelliti et al., 2005). In this study, E _f was set to −20 mV, and G _f was chosen to be 0.2 and 2 nS, to investigate the effects of various parameters in the cellular fibroblast model on the output of the simulations.

2.1.2.  Myocyte-fibroblast mechanical modeling at the cellular level: the Ca²⁺-force relation

 To simulate the mechanics of the myocyte-fibroblast model, the Ca²⁺-force relation in myocytes and fibroblasts has been considered. In this study, the Rice mechanical model of a cardiac myofilament (Rice et al., 1999; Iribe et al., 2006) was applied for the myocyte. It was coupled to the modified ten Tusscher model by using the Ca²⁺ concentration as the input to the Rice model.

 To the best of our knowledge, there is no well-defined mathematical model to describe the tension in fibroblasts. The main reasons are that: (1) There is lack of experimental data due to the small size of individual fibroblasts (Brown et al., 1998); (2) Ca²⁺ current, as a control factor of active tension, has not been identified for membrane ionic currents in fibroblasts (Chilton et al., 2005; Shibukawa et al., 2005; MacCannell et al., 2007; Maleckar et al., 2009). The active force of fibroblasts was therefore ignored at the cellular level in this study.

2.1.3.  Tissue modeling

 With the above mentioned cellular electromechanical models for myocytes and fibroblasts, both AP and T _a can be computed and then incorporated into the following tissue models.

 At the tissue level, the fibroblasts and myocytes are assembled into a 2D tissue sheet and the electrophysiological and mechanical coupling is calculated to simulate the excitation conduction and deformation in the heart.

 The equation of myocyte-fibroblast electrophysiological modeling is expressed as (Xie et al., 2009a) (5) , where C is the capacitance of the myocyte (C _m) or fibroblast (C _f); V is the membrane potential of the myocyte or fibroblast; I _ion is the current of myofilament (I _m) or fibroblast (I _fb); n is the number of the coupled neighboring cells (either myocyte or fibroblast cells); and, is the gap junction conductance between a cell (either a myocyte or a fibroblast cell) and its kth neighbor (either a myocyte or a fibroblast cell). According to the experimental reports, C _f ranges from 6.3 to 75 pF (Vasquez et al., 2004; MacCannell et al., 2007). The G _gap between the myocytes ranges from 262 to 937 nS for end-to-end coupling (Yao et al., 2003), and the G _gap between a myocyte and a fibroblast or between two fibroblasts (G _{gap_j}) ranges from 0.3 to 8 nS in cultured cells (Rook et al., 1992), whereas 0–100 nS were used in other modeling studies (MacCannell et al., 2007; Jacquemet and Henriquez, 2008; Sachse et al., 2008). In this study, G _gap and C _f were set to 500 nS and 25 pF, respectively. G _{gap_j} was chosen as 4 and 20 nS. From Eq. (5), the excitation conduction of the 2D cardiac tissue can be realized.

 For tissue mechanics, the mechanical model of cardiac tissue should contain the mechanics of normal myocyte tissue and fibroblast proliferous tissue. The mechanical model proposed by Nash and Panfilov (2004) was used to model the mechanics of both tissues, and includes the stress equilibrium, the finite element approximations, the constitutive law, and the uniform isotonic boundary loads. The main equations are expressed as follows: (6) , (7) , (8) , (9) , where X_M is the material coordinate; T^MN is the second Piola-Kirchhoff stress tensor of the finite elasticity theory; is the deformation gradient tensor; W is the scalar strain energy density function; E_MN is the Green’s strain component; C^MN is the component of the contra-variant metric tensor; Ψ_n is the geometric interpolation function; V ₀ is the undeformed volume; S ₂ is the portion of the boundary subject to external tractions; f_n ^j is the external nodal traction force; n_j is the coordinate of a vector normal to the boundary; and, p is the uniform load.

 Because fibroblast proliferation occurs after a myocardial infarction and cardiac fibroblast proliferous remodeling leads to a progressive increase in ventricular passive stiffness (Conrad et al., 1995; Biernacka and Frangogiannis, 2011), the common approach for modeling ischemic myocardium is to change the material descriptions of the passive myocardium. In this study, the material constants in the passive material properties of fibroblast proliferous tissue were re-evaluated to show the higher elastic stiffness of the myocardial scar tissue. Compared to the normal tissue with material constants c ₁=2 kPa and c ₂=6 kPa (Nash and Panfilov, 2004), c ₁ and c ₂ in the fibroblast proliferous tissue were set to 4 and 12 kPa, respectively. In this way, the modulus of the proliferous tissue rose to double that of the normal tissue, within the range of the experimental results (Berry et al., 2006). In addition, the central and edge nodes were completely fixed in space to prevent rigid-body translations and rotations. After the calculation of the mechanical equations, the deformation of the 2D tissue can be calculated.

 The main framework of the proposed coupling model has been described in detail. It integrates two kinds of models (electrophysiological and mechanical models), two kinds of cells (myocytes and fibroblasts), and two spatial levels (cellular and tissue levels). The proposed model includes the sum of the ordinary differential equations (ODEs) and partial differential equations (PDEs), which requires a large computational effort. T _a and the electrophysiological model of a cardiac myocyte and a fibroblast were represented by ODEs, and the reaction diffusion of cardiac conduction and the stress equilibrium of cardiac mechanics were represented by PDEs. In the following section, the computational procedures are given in detail.

2.2.  Computational procedures

 In the 2D coupling model, we used tissue sizes of 225×225 grid points to study the CV of plane wave propagation, depolarization and repolarization times, and the APD of a point stimulus (Cherry and Fenton, 2011), and 603×603 grid points to study spiral wave rotation. Owing to the mixed myocyte-fibroblast coupling tissue, we used different diffusion tensors based on different situations of coupled neighboring cells: the ‘diffusion’ coefficient D=1.24 cm²/s to simulate myocyte-myocyte coupling and D=0.94 cm²/s (G _{gap_j}=20 nS) or 0.81 cm²/s (G _{gap_j}=4 nS) to simulate myocyte-fibroblast or fibroblast-fibroblast coupling (ten Tusscher and Panfilov, 2005). For the finite difference method (FDM), the spatial resolution (Δx) was 0.1 or 0.2 mm, and the temporal resolution (Δt) was 0.005 or 0.010 ms. No-flux boundary conditions were used. For the finite element method (FEM), we used a mesh with 14×14 elements containing 16×16 grid points in plane wave and central point stimulus simulations, and 43×43 grid points in the re-entry simulation. The temporal resolution for the mechanical model was 2 ms. Strain maps of the central stimulus and re-entry, and temporal traces of strain in different points of the tissue were also investigated.

 The system was solved numerically using custom software written in Fortran language. At the cellular level, all state variables were updated by the forward Eular method. The FDM was applied to solve the reaction-diffusion equation. Following each time integration step, all parameters of cells were updated. The active stress was then interpolated at the four-node rectangle isoparametric element Gaussian points. The stresses of these active Gaussian points served as the inputs to the governing equations of the tissue mechanics model. The stress equilibrium equation was solved by a nonlinear least square iteration method with different material constants from different tissues.

 In simulations, the electromechanical properties of a single human myocardial cell and a fibroblast were first computed using the proposed integrative model, and were then compared to the original characteristics of the corresponding electrophysiological and mechanic models—ten Tusscher and the Rice models. Second, the influence of diffuse fibroblast proliferation on plane wave propagation in 2D tissue was investigated. Third, the electromechanical coupling of a central point stimulus with two spatial resolutions, two temporal resolutions, two gap junction conductances, and two fibroblast membrane conductances was simulated and compared, to investigate the effect of various parameters on depolarization and repolarization times, APD, strain maps, and temporal traces of strain in different points of the tissue. Finally, the electromechanical coupling of the re-entry and strain maps was simulated. The spatial configuration of normal tissue and fibroblast proliferous tissue is illustrated in Fig. 2.

3.  Results

3.1.  Electromechanical coupling at the cellular level

 The myocardial-fibroblastic electromechanical coupling at the cellular level was first simulated and then compared to the original ten Tusscher model. The AP, Ca_i, and T _a curves of the myocyte generated by the proposed coupled model after 10 stimuli are shown in Fig. 3.




 The AP and Ca_i computed from the original ten Tusscher model were also compared. The main parameters of both models are shown in Table 1. Note that the maximum Ca_i in the coupled model was lower than that in the ten Tusscher model (Fig. 3). This was caused by the modification of the cytosolic concentration of Ca²⁺ (Eqs. (1)–(3)). The concentrations of the Ca²⁺ bound to calmodulin and troponin were removed from the cytosolic Ca²⁺ concentration and hence reduced Ca_i, compared to the value obtained experimentally by Beuckelmann et al. (1992). The calcium dynamics affected Ca_i, which further affected the APD of the coupled model. In Fig. 3c, T _a has a maximum force of 0.0205 N/mm² (20.5 kPa), which is consistent with the general range (10–75 kPa), depending on the sarcomere length (Niederer et al., 2006). Collectively, the proposed model can well reproduce the electrical AP of the original ten Tusscher model, validating the proposed coupling model to a certain degree.




 As given in Eq. (5), the fibroblast-myocyte coupling at the tissue level can be calculated with various numbers of fibroblasts coupled to a myocyte. Here, we coupled one ventricular myocyte to three passive fibroblasts. The AP and T _a obtained after 10 stimuli are shown in Fig. 4.




 Compared to the AP of the uncoupled myocyte, the peak of the myocyte AP of the coupled model decreased from 36.3 to 26.5 mV, with the coupled ventricular APD increasing from 289 to 335 ms (Fig. 4a). The resting potential depolarized slowly within 20 ms and increased by 4.9% at 20 ms when compared with the uncoupled control. The peak of T _a was slightly decreased by 3.9% with the coupled fibroblasts (Fig. 4b). These simulation results demonstrate that fibroblasts do affect cardiac electromechanics at the cellular level.

 Fig. 5 presents the APD restitution curve at 90% repolarization and the half-width (HD) of the averaged Ca_i transient from a control myocyte vs. a myocyte coupled to three passive fibroblasts. The pacing rate ranges from 0.5 to 2.5 Hz. The longer APD and HD of fibroblast-myocyte coupling can be observed.

3.2.  Planar conduction velocity

 The influence of diffuse fibroblast proliferation on plane wave propagation was investigated in the first set of simulations. The Δx was set to 0.1 mm and Δt was set to 0.01 ms. In Figs. 6a and 6b, waves were initiated at the lower border of the cardiac tissue with no or 30% diffuse fibroblast proliferous tissue, respectively. Snapshots were taken at the same time instant. Fig. 6c shows the planar CV as a function of the percentage of fibroblast proliferous tissue. With different values of G _{gap_j} and G _f, four different cases of CV are illustrated. The lower G _{gap_j} resulted in an earlier propagation block. When G _{gap_j}=4 nS, a propagation block occurred when fibroblast proliferation reached more than 40%, compared with 45% when the G _{gap_j}=20 nS. The value of G _f also had a slight impact on the CV.

3.3.  Two-dimensional simulation with a central stimulus

 In this set of simulations, electromechanical activities and the change of strain caused by a periodic point stimulation at the center of the tissue were investigated. Central stimulus sites were chosen to be the smallest possible square regions (7×7 computational nodes), which made it possible to produce a propagating wave. Ten periodic central stimuli with a time interval of 380 ms were applied at the center of the 2D tissue to ensure a stable excitation and contraction; the 11st activation was recorded. Other parameters were listed as follows: Δx=0.1 mm, Δt=0.005 ms, G _{gap_j}=4 nS, and G _f=0.2 nS.

 Figs. 7a and 7c show snapshots of the electromechanical activities in the 2D planes of cardiac tissue with no or 30% diffuse fibroblast proliferous tissue, respectively, in three areas as shown in Fig. 2. Figs. 7b and 7d show the corresponding deformations. The emergence of T _a lagged behind the excitation of the myocyte (Figs. 7a and 7c). Hence, the deformations introduced by the mechanical contraction in the 2D tissue also slightly lagged behind the radically-propagating waves of the excitation. At 3810 and 4200 ms, the fibroblast proliferous areas slowed the propagating waves. In addition, the meshes remained slightly relaxed in the center region when the excitation waves spread outward. This was caused by the last cycle of the mechanical contraction and the inertia of the elastic tissue. At 3910 and 4380 ms, the contractions were very clear. Because of the higher elastic stiffness of the myocardial scar tissue, the deformations of the fibroblast proliferous tissue were significantly lower than that of the non-proliferous tissue, about 0.05 vs. 0.1. Throughout the whole process, the deformations in the non-proliferous tissue elicited by the propagating waves of excitation were symmetric about the central axes (Fig. 7b), while being non-symmetric in fibroblast proliferous tissue (Fig. 7d). Note that for both tissues, deformations were highest in the central regions and gradually decreased from the inside to the outside. Moreover, comparing the peak values of deformations in these two tissues, we can see that as fibroblast proliferous tissue was added, deformations in the non-proliferous tissue increased, especially in the central area, from 0.18 to 0.23.

3.4.  Spatial profiles in square geometries with a central stimulus

 To explore the possible effects of varying parameters on simulation results, different combinations of two spatial resolutions (Δx=0.1 or 0.2 mm), two temporal resolutions (Δt=0.005 or 0.01 ms), two gap junction conductances (G _{gap_j}=4 or 20 nS), and two fibroblast membrane conductances (G _f=0.2 or 2 nS) were studied. In this study, two groups with different Δt are presented. Fig. 8 shows the depolarization, repolarization, APD, and histograms of APD dispersion profiles from the normal and fibroblast proliferous tissues. Dispersions of the depolarization, repolarization, and APD in Fig. 8 were given in Table 2. In this set of simulations, the mechanoelectric feedback was taken into account. In other words, excitation conduction was related to tissue deformation.







 From Fig. 8 and Table 2, it can be seen that: (1) All the dispersions in the depolarization, repolarization, and APD increased after a fibroblast proliferation presented in the tissue. (2) For each Δt, increases in dispersions of the depolarization in the two tissues were 8 and 10 ms, these increases were greater than those of the repolarization (5 ms) and APD (2 ms). (3) With increasing Δt, dispersions in the depolarization rose from 32 to 45 ms in the normal tissue, and from 40 to 55 ms in the fibroblast proliferous tissue, increases of 13 and 15 ms, respectively. Similarly, dispersions in the repolarization in the two tissues increased by 15 ms. Dispersions in APD were not sensitive to Δt. In addition, the depolarization and repolarization time in the three proliferous regions lagged behind those in the corresponding locations in the normal tissue with no fibroblast proliferation (Fig. 8). This further indicated that fibroblast proliferation slowed down the depolarization and repolarization, because of the small gap junction conductance and fibroblast membrane conductance.

 The histograms of APD dispersions shown in Fig. 8 were measured radially every 5° around circular wavefronts. For each Δt, APD dispersions in the normal tissue changed uniformly from 0° to 360°. However, in the fibroblast proliferous tissue, APD dispersions along directions from the center to proliferous regions increased significantly.

3.5.  Temporal trace of strain at different points and histograms of Point 4 with a central stimulus

 Fig. 9a shows temporal traces of the strain at four points in the normal tissue with Δt=0.005 ms. The closer they were to the central site of the 2D tissue, the larger the strain they sustained. Point 4 was the closest to the central site and therefore showed higher strain than the other three points.




 Fig. 9b shows histograms of the maximum strain and HD at Point 4 with different spatial and temporal resolutions, gap junction conductance, and fibroblast membrane conductance. First, for any Δt and Δx, fibroblast proliferation led to an increase in all the maximums and HDs. Second, for the same Δt, an increase in Δx led to decreases in strain peaks and increases of HDs. For the same Δx, the increase of Δt led to the same trend in strain peaks and HDs. Third, changes in G _{gap_j} and G _f had little impact on the strain. The determining factor was the modulus of the fibroblast proliferous tissue.

3.6.  Two-dimensional simulation of spiral waves

 In this set of simulations, electromechanical activities and a change of strain with re-entrant excitation waves were performed. Re-entrant excitation waves were initiated using a standard S1-S2 protocol (ten Tusscher et al., 2004). First, a single S1 stimulus was applied to the left side of the tissue at 20 ms after the initial condition, to produce a planar wave that propagated rightward. Then, a second S2 stimulus was applied to the lower-left quarter portion of the tissue when the left-half portion of the tissue returned to its resting state. The S1 and S2 stimuli had durations of 1 and 5 ms, respectively, but both had twice the amplitude of the diastolic excitation threshold (ten Tusscher et al., 2004). Other parameters were: Δx=0.2 mm, Δt=0.01 ms, G _{gap_j}=4 nS, and G _f=0.2 nS.

 Figs. 10a and 10c show snapshots of electromechanical activities in 2D planes of cardiac tissue with no fibroblast proliferous tissue and with 30% of diffuse fibroblast proliferous tissue in three areas (Fig. 2). Figs. 10b and 10d show the corresponding deformations. In Fig. 10c, the tip of the spiral wave was trapped to rotate around the border of the fibroblast proliferous zone at the middle left, rather than at the center, as shown in Fig. 10a. In the normal tissue with no fibroblast proliferation, the rotor of the spiral wave had been drifting to the left, while in the fibroblast proliferous tissue, it moved around the fibroblast proliferous region with a smooth rotary action. This indicated that fibroblast proliferation can sustain stable spiral wave rotation.




 In Fig. 10b, the deformation in the rotor was larger than in the other regions. As time went on, the deformation gradually increased. In the fibroblast proliferous tissue, larger deformations did not follow the tip of the spiral wave, but focused on areas among the fibroblast proliferous regions (Fig. 10d). Owing to the slower propagation of the excitation wave in the fibroblast proliferous tissue, these regions had long been in the excited state, which resulted in a sustained contraction in the three fibroblast proliferous areas.

4.  Discussion

 In this study, we developed a strongly coupled myocardial-fibroblastic electromechanical model to explore the functional effects of fibroblast proliferation, including myocyte-fibroblast couplings, on the modulating cardiac AP morphology at the cellular level, and cardiac excitation wave conduction and contraction at the tissue level. The influence of diffuse fibroblast proliferation on plane wave propagation was studied. Parameters were varied, such as the spatial and temporal resolutions, the conductance of the fibroblast and the gap junctional coupling, to investigate their effects on the spatial distribution of depolarization, repolarization, and APD. Strain maps and temporal traces of strain at different points in both normal and fibroblast proliferous tissues were also illustrated. Electromechanical models of a central point stimulus and re-entry were simulated. The major findings of this study were: at the cellular level, the myocyte-fibroblast coupling prolonged the APD; at the tissue level, the CV decreased as fibroblast proliferation increased and a conduction block occurred at 40%–45% fibroblast proliferation. For the central stimulus, owing to the higher elastic stiffness of the myocardial scar tissue, the deformations of fibroblast proliferous tissue were significantly lower than that of non-proliferous tissue. The dispersions in the depolarization, repolarization, and APD all increased after a fibroblast proliferation presented in the tissue. For re-entry, the fibroblast proliferation sustained spiral waves and kept the contraction in the fibroblast proliferous regions.

4.1.  Reasonableness of the model assumptions

 Anatomic analysis indicates that the heart is composed of myocytes, vasculature cells, and connective tissue cells (Manabe et al., 2002). However, previous mathematical models of cardiac strong/weak electromechanical couplings had not considered fibroblasts (Xia et al., 2005; Kerckhoffs et al., 2006; Trayanova et al., 2011). In this study, we explored the possible influence of fibroblasts on cardiac electrophysiology and mechanics by adding fibroblasts to a 2D electromechanical cardiac tissue model. Diffuse fibroblast proliferation was modeled by the presence of size 1×1 grid points that were randomly distributed in three areas of the 2D tissue. The three regions of fibroblast proliferous tissue were set to be flexible structures, but with a large elastic modulus. This assumption was reasonable: cardiac fibroblast proliferous remodeling leads to a progressive increase in passive ventricular stiffness (Conrad et al., 1995; Berry et al., 2006; Biernacka and Frangogiannis, 2011), and the flexibility in these tissues is very poor, especially in cardiac fibroblast proliferous tissue (Brown et al., 1998). Therefore, local strains should be relatively small within the fibroblast proliferous regions. However, they can still be displaced by forces acting from the surrounding myocardium.

4.2.  Effects of fibroblast coupling on human ventricular myocyte AP waveform and T _a

 For fibroblasts, previous studies using myocyte-fibroblast co-culture models have shown that the fibroblasts depolarized electrotonically coupled myocytes (Miragoli et al., 2006; MacCannell et al., 2007; Jacquemet and Henriquez, 2008; Maleckar et al., 2009). Our simulations also showed a depolarizing effect of coupled fibroblasts on the resting membrane potential of ventricular myocytes.

 Coupling the human ventricular myocyte with fibroblasts in this study also modulated APD. In previous studies, both “passive” and “active” fibroblasts resulted in diverse effects on APD. Xie et al. (2009a) used a passive model and showed that myocyte-fibroblast coupling prolonged or shortened APD based on the properties of G _f and E _f. MacCannell et al. (2007) and Jacquemet and Henriquez (2008) used the same active model but found the opposite tendencies in APD, depending on the resting potential of the fibroblasts. In our simulations, the myocyte-fibroblast coupling prolonged the APD by 46 ms, which accorded with previous modeling studies (Jacquemet and Henriquez, 2008; Xie et al., 2009a).

 In addition to the effect on human ventricular myocyte electrophysiology, coupled fibroblasts had implications for changes in T _a. Kerckhoffs et al. (2010) showed that in the failing heart model, inotropy was decreased by reducing peak fiber T _a by 27%. In our study, the peak of T _a decreased by 3.9% with the coupled fibroblasts. The reason was that fibroblasts coupling reduced the Ca²⁺ concentration and caused a reduction in T _a. Due to the strong correlation between fibroblast proliferation and heart failure (John et al., 2004; Assomull et al., 2006; Pellman et al., 2010), it can be expected that once the number of fibroblasts increases, a significant reduction would occur in T _a.

4.3.  Effects of fibroblast proliferation on CV

 In an experimental study, Miragoli et al. (2006) found a biphasic effect on CV when endogeneous fibroblasts proliferated with myocytes. In a modeling study, Xie et al. (2009a) showed that fibroblasts, whether coupled to myocytes or not, slowed conduction by creating zigzag conduction pathways. Our simulation showed that G _{gap_j} and G _f had significant effects on wave propagation. A smaller G _{gap_j} caused a conduction block at a lower percentage of fibroblast proliferation. A smaller G _f made the CV drop relatively quickly. Qualitatively, the results of CV in this study were substantially different from previous findings using the FitzHugh-Nagumo model and the ten Tusscher human ventricular tissue model (ten Tusscher and Panfilov, 2003; 2007). In these models, the CV decreased more strongly as a function of the percentage of fibroblast proliferation (25% and 21% CV decreases, respectively, for 20% fibroblast proliferation) than in our simulations (13.3%–16.7% CV decrease for 20% fibroblast proliferation, based on different values of G _{gap_j} and G _f). This was because the diffuse fibroblast proliferous tissue was modeled by unexcitable obstacles in previous models, while in our study, the tissue could still conduct excitement with a small gap junctional conductance.

4.4.  Effects of fibroblast proliferation at the tissue level

 For central stimulus simulations, the effects of fibroblasts on cardiac excitation conduction and contraction were investigated in several ways. Our simulations indicated that changes in G _{gap_j} and G _f were not the determinate factors of the strain, but the modulus of the fibroblast proliferous tissue. Antoni et al. (2010) showed that strain in normal regions was (−15.3±4.5)%, and (−10.6±5.3)% in the infarct zone. In our calculations, peak deformation in normal tissue was 18.0%, which was in the range, while in fibroblast proliferous tissue, it was up to 23.0%, which was slightly higher than the data. Based on these results, we suspect that once severe fibroblast proliferation occurred, strains decreased, apparently because of the large stiffness in these areas or because of the similarity to rigid motion. It would be interesting to check this prediction against available clinical data. Cherry and Fenton (2011) discussed the importance of tissue size in determining the spatial distribution of APD in normal homogeneous tissue and showed that for the small and large squares, the magnitude of the dispersion of depolarization was the largest, while the total dispersion in APD was nearly identical. Our simulations verified these results and also proved that increases in depolarization, repolarization, and APD dispersions accompanied the occurrence of fibroblast proliferation.

 For spiral wave simulations, once fibroblast proliferation presented in the 2D tissue, the spiral wave did not propagate along its inherent directions and had to bypass fibroblast proliferous tissues. In this way, the fibroblast proliferation changed the form of the re-entry. In addition, because of the presence of the fibroblasts, the CV slowed (Miragoli et al., 2006; Zlochiver et al., 2008), which contributed to stabilization of the spiral wave. This further illustrates that fibroblast proliferation aggravated the cardiac condition. For deformations, owing to the sustained excitement in the core of the spiral wave, active stresses in the rotor remained at a high level, which made the strain become larger than that in other regions. In fibroblast proliferous tissue, because the propagation of the excitation wave in the fibroblast proliferous regions was slower than that in normal tissue, these regions had long been in the excited state, which resulted in a sustained contraction in the three fibroblast proliferous areas.

4.5.  Limitations

 Several limitations of this study should be mentioned. First, mechanosensitive currents have not been incorporated into the coupled model. Second, due to the lack of physiological data from the same species for the coupled model, only the existing experimental data relevant to each individual model were used to validate our simulation results. Third, the mechanisms of the fibroblast mechanics should be further studied. Previous studies of fibroblast mechanics were based mainly on experiments. Once a detailed mechanical model of the fibroblast is developed, it can be coupled to the myocyte mechanics model to simulate more precisely the mechanisms of the myocardial-fibroblastic electromechanical coupling.

5.  Conclusions

 A strongly coupled myocardial-fibroblastic electromechanical model has been developed by integrating a fibroblast model with a cardiac electrophysiological and mechanical model, and the effects of fibroblasts on cardiac excitation conduction and contraction have been investigated. The simulation results showed that the fibroblasts slowed wave propagation, induced a conduction block, decreased strains in fibroblast proliferous tissue, increased the dispersions in depolarization, repolarization and APD, changed wave forms, and sustained re-entry. These identified effects demonstrate that fibroblasts are important in electromechanical coupling and should be considered in future cardiac electromechanical modeling.

---

* Project supported by the National Basic Research Program (973) of China (No. 2007CB512100) and the National Natural Science Foundation of China (Nos. 81171421 and 61101046)Compliance with ethics guidelines He-qing ZHAN, Ling XIA, Guo-fa SHOU, Yun-liang ZANG, Feng LIU, and Stuart CROZIER declare that they have no conflict of interest.