Abstract
The purpose of the paper is to derive a formula to describe the quantitative relationship among the number of the pacing cells required (NPR), the dimension , and the diffusion coefficient (electrical coupling or gap junction ). The relationship between NPR and has been investigated in different dimensions, respectively. That is, for each fixed , there is a formula to describe the relationship between NPR and ; and three formulas are required for the three dimensions. However, there is not a universal expression to describe the relationship among NPR, , and together. In the manuscript, surveying and investigating the basic law among the existed data, we speculate the preliminary formula of the relationship among the NPR, , and ; and then, employing the cftool in MATLAB, the explicit formulas are derived for different cases. In addition, the goodness of fit () is computed to evaluate the fitting of the formulas. Moreover, the 1D and 2D ventricular tissue models containing biological pacemakers are developed to derive more data to validate the formula. The results suggest that the relationship among the NPR, , and the () could be described by a universal formula, where the NPR scales with the (the dimension) power of the product of the square root of () and a constant which is dependent on the strength of the pacing cells and so on.
1. Introduction
In the normal heart, the electrical pulses are initiated in the genuine pacemaker–sinoatrial node (SAN), which generates the excitation automatically [1, 2]. However, the dysfunction of the SAN would lead to a variety of manifestations, including fibrillation, arrhythmias, and heart failure [3–5]. To treat these diseases, the best way at present is to implant electronic devices [6]. There are more than 200,000 patients suffering electrical pacemaker implantation every year in USA alone [7]. Nevertheless, the limitations of the devices could not be ignored, including lead malfunction, infection, short battery lifespan, and thrombosis and so on [8–12].
As a consequence, the biological pacemaker has been attracting the attention of the researchers to overcome the disadvantages of the electronic devices [13–16]. One of the popular strategies is to create biological pacemakers in the ventricle [17–19], because the thicker ventricular wall is more conductive to the biological experiments.
For a successful pacemaker, an important feature is the source-to-sink match. The pacemaker acts as the source to drive the adjacent resting cardiac tissue which is considered as the sink. The sink is at a resting state until it is stimulated by the pulses from the source. That is, the source must be strong enough to drive the sink to break the threshold to depolarize. The depolarization of the sink would fail if the source-to-sink mismatch is too large, leading to the failure of the excitation propagation. The topic has been investigated extensively in the SAN [20, 21]. The research demonstrates that the gap junction (electrical coupling) plays an important role in the source-to-sink mismatch [22–24]. The coupling conductance is much weaker in the SAN than that in other cardiac tissues [25, 26]. The poor coupling reduces the suppression from the adjacent hyperpolarized tissue, shielding the depolarization of the pacing tissue to overcome the source-to-sink mismatch [27].
The ability of the source (the pacing tissue) to drive the sink (the cardiac tissue) depends on the number of the pacing cells, the gap junction conductance, and the strength of the pacing cells. For a successful pacing and driving, what is the relationship between the NPR and the gap junction conductance? Tveito et al. showed that the NPR was proportional to the square root of the gap junction conductance in the1D strand and to the gap junction conductance in the 2D tissue [28, 29]. Xie et al. validated the consistent conclusion for the propagation of the action potential triggered by early (EAD) and delayed afterdepolarizations (DAD), and they further concluded that the NPR scaled with the 1.5 power of the gap junction conductance in the 3D volume [30].
The analysis of the relationship between the NPR and the gap junction conductance to overcome the source-to-sink mismatch has been performed in different dimensions, respectively. That is, for each fixed dimension , there is a formula to describe the relationship between NPR and the gap junction , and three formulas are required for all the three dimensions. However, to our knowledge, there is not a universal expression to show the precise relationship among the NPR, gap junction , and dimension together. In the manuscript, using the computational model to obtain the simulation data, and together with the data from reference [30], we set out to further probe the topic and try to address a universal formula to show how many pacing cells in different dimensions are required to overcome the source-to-sink mismatch to make the pulses propagate from the pacing tissue with variable gap junction conductance into the surrounding regions.
Firstly, surveying and analyzing the existing data [30] on the whole, we speculate and derive a basic format of the formula describing the relationship among the NPR, , and . And then, the accurate formulas were calculated, validated, and evaluated for different cases. Thirdly, based on the TNNP06 ventricular single-cell model [31], the pacing cell model is addressed by depressing to 0.05 nS/pF. Next, using the reaction diffusion equation, the 1D and 2D models are developed by coupling the pacing cells and ventricular myocytes together. From the models, the NPRs are obtained corresponding to variable gap junction conductance in different dimensions. Finally, the simulation data are utilized to verify the formula derived previously.
2. Data and Methods
In this section, we introduce the models of the single cell and the tissue from which the simulation data are derived; and then, the related data from other reference are listed.
2.1. The Model of the Single Cell
For the single cell, we employ the well-established TNNP 2006 model of the human ventricular cells [31], which is shown in Equation (1). where
In Equation (1), is the transmembrane potential; is the sum of all the transmembrane ion currents; is the external stimulus current and is set to 0 in the study; is membrane capacitance per unit surface area; is the inward-rectifier K+ current; ,, and are the outward slow, rapid, and transient and rectifier potassium currents, respectively.
In the simulation, is depressed to 0.05 nS/pF to obtain the robust pacing cells.
2.2. The Models of the 1D and 2D Tissues
The settings of the 1D and 2D tissue are as shown in Figures 1(a) and 1(b), respectively. The central part, consisting of pacing cells, is the pacemaker, and the surrounding blue region is the ventricular tissue which is made up of endomyocardial myocytes.
(a)
(b)
The reaction-diffusion equation [31] is applied to describe the propagation of the electrical excitation wave in the 1D and 2D ventricular tissues, which is shown in Equation (3): where is the diffusion coefficient with the normal value 15.4 mm2/s, which is corresponding to the electrical coupling (gap junction) in the tissue and describes the conductivity of the excitation; is the Laplace operator; and the other parameters are the same as those in Equation (1). In the simulation, the time step is set to 0.02 ms, and the space step is 0.15 mm.
2.3. The Data from the Existing Reference
Xie et al. investigated how the source overcame the mismatch to trigger the successful pulse propagation in different conditions [30]. The numbers of contiguous automatic cells required are listed in Tables 1 and 2 for various sources with diverse gap junctions in different dimensions.
In Tables 1 and 2, the abbreviations are listed in Table 3.
“EAD/normal,” “EAD/RRR,” “DAD/normal,” and “DAD/HF” are four pacemaker-tissue designs (cases).
2.4. The MATLAB Tool
The cftool in MATLAB is adopted for the fitting in the manuscript.
3. Results and Discussion
In this section, our aim is to derive a formula to fit the relationship among the NPR, the electrical coupling (gap junction ), and the dimension And then, the accuracy of the formula is evaluated.
3.1. The Fitting Formula
From the observation and computation, we detect that in Tables 1 and 2, for each column of data, is approximately a geometric sequence. The details are shown in Tables 4 and 5.
For all the 8 columns of data in Tables 4 and 5, is close to , which infers that for a fixed , might be geometric sequences corresponding to , and is the common ratio . As a consequence, the expression of of can be described as where and are functions of , and is the dimension.
What is more, the average values of for each in different cases are computed according to the values of in Tables 4 and 5, respectively, which are listed in Table 6.
From Table 6, for all the cases “EAD/normal,” “EAD/RRR,” “DAD/normal,” and “DAD/HF,” it is found that is close to , that is, is proportional to . As a consequence, we set , where is a constant.
Therefore, the expression of about can be rewritten as Equation (5).
According to [28–30], for a fixed , is proportional to the power of the square root of the gap junction (coupling conductance) in the dimensional tissue, and therefore, could also be expressed as follows: where is a function of , and is the gap junction (coupling) conductance.
Considering Equations (5) and (6) together, for and , there is
Therefore, must be a constant, set as . And the formula of could be described as Equation (8). where and are the undetermined coefficients.
In the previous studies [28–30], only the relationship between NPR and for each fixed dimension is found. As a consequence, there are different formulas for different dimensions. That is, there are three formulas for all the three dimensions and three parameters should be determined. In the study, the relationship between NPR and for each fixed is also discussed. Based on the two relationships, Equation (8) is obtained. And only one formula is required for all the three dimensions and only two parameters are needed.
Using the cftool in MATLAB, we set out to fit the formula for the case “EAD/normal” according to the first two columns of data in Table 1. Knowing the format of the formula in Equation (8) in advance, we derive the formula for the case “EAD/normal”:
The computing for the case “EAD/normal” according to Equation (9) are shown in Table 7, compared with the original values.
The comparison in Table 7 illustrates that all the computing results are close to the corresponding original values; that is, the function fits the original data well.
The graph of the fitting function is presented in Figure 2 together with the 6 original points in the first two columns of the data in Table 1.
From Figure 2, it could be seen intuitively that the small circles (the original values) cling to the surface of the function of , which means that formula (9) is a good fitting.
The analogous operations are done for the cases “EAD/RRR,” “DAD/normal,” and “DAD/HF.” The formula for case “EAD/RRR” is shown in Equation (10).
And the graph of function (10) is presented in Figure 3 together with the original values of the case “EAD/RRR” in the second two columns of data in Table 1.
The formula for case “DAD/normal” is shown in Equation (11).
And the graph of function (11) is presented in Figure 4 together with the original values of case “DAD/normal” in the first two columns of data in Table 2.
The formula for case “DAD/HF” is shown in Equation (12).
And the graph of function (12) is presented in Figure 5 together with the original values of case “DAD/HF” in the second two columns of data in Table 2.
3.2. The Evaluation of the Formula
In the previous section, we show the intuitive fitting results in Figures 2–5, which demonstrate that the formula in Equation (8) fits all the 4 cases accurately.
For a fixed , the formula is an exponential function of , and it could be transformed logarithmically to derive a linear function. And then, the goodness of fit (), which is an excellent tool to evaluate the linear function, is calculated to evaluate the fitting.
Taking log of both sides of Equation (8), we can obtain
And the linear expression for the case “EAD/normal” is
The graphs of function (14) for and are shown in Figure 6. And the 6 logarithmic original in Table 1 for case “EAD/normal” are drawn in pink for and in red for , respectively.
Figure 6 illustrates that the formula fits the original values well intuitively. Furthermore, is computed to evaluate the fitting quantitatively.
The closer is to 1, the better the fitting is. And values are 0.9964 and 1.0000 for and for , respectively. That is, the fitting for case “EAD/normal” is accurate and acceptable.
The linear transformation for cases “EAD/RRR,” “DAD/normal,” and “DAD/HF” is expressed in Equations (15), (16), and (17), and the corresponding graphs are shown in Figures 7–9.
The goodness of fit () for all the four cases are listed in Table 8.
Intuitively, all the lines in Figures 7–9 fit the corresponding original values closely. Moreover, all the values are no less than 0.99, inferring that the fitting are accurate.
3.3. The Verification of the Formula
In this section, using the models in Section 2.1 and Section 2.2, we simulate the biological pacemakers with different couplings to derive the corresponding number of the pacing cells required. And then, we set out to verify whether formula (8) is suitable for the case.
The structures of the 1D and 2D tissues are as shown in Figure 1. The pacemakers are in the central region, composed of pacing cells which are transformed from the ventricular myocytes by depressing to 0.05 nS/pF. The radius of the pacemaker is increased progressively until the pulses generated automatically from the pacemaker propagate successfully to the surrounding ventricular tissue.
The NPRs corresponding to the different diffusion coefficient are shown in Table 9.
Next, the relationship between gap junction and diffusion coefficient is introduced first. The diffusion throughout a coupled tissue depends on the presence of the gap junction . And is also a tensor describing the conductivity of the tissue.
In the study, Equation (3) is adopted from reference [31] to describe the propagation of the electrical excitation wave in the cardiac tissue, which could be discretized as Equation (18). where dh is the space step and is set to 0.15 mm in the simulation. What is more, is set to 0 in the paper, because the electrical excitation is generated from the pacemaker and the external stimulus is not required. As a consequence, the final discretized reaction diffusion equation adopted in the paper is Equation (19).
On the other hand, the data in Tables 1 and 2 are cited from reference [30], which are derived by Equation (20) describing the propagation of the electrical excitation.
Comparing Equation (19) with Equation (20), it is found that is corresponding to . That is, except for a coefficient, and are essentially the same. As a consequence, the data about NPR, , and will be suitable to evaluate the formula (8).
Based on formula (8), we perform the analogous analysis using the cftool in MATLAB and derive the quantitative relationship among the NPR (), the diffusion coefficient , and the dimension , described in Equation (21).
The graph of the function is shown in Figure 10 together with the 12 original values in Table 9 to intuitively evaluate the fitting of the function (21).
Figure 10 illustrates that the original values are close to the surface of the function, which suggests that the function fits the data well.
According to formula (21), when , the equation could be rewritten as
Based on formula (22), the graph is shown in Figure 11 together with the corresponding original values in Table 9.
And, when ,
The graph of function (23) is shown in Figure 12 together with the corresponding original values in Table 9.
The goodness of fit () for formula (22) and formula (23) are 0.9998 and 0.9999, respectively, which demonstrates the fitting is accurate.
In summary, based on the relationship between NPR and in references [28–30] and the relationship between NPR and found in the study, together with the data from reference [30] and the simulation in the paper, the relationship among the NPR (), the dimension , and the gap junction (diffusion coefficient ) is derived and evaluated, which satisfies the two-variable function in Equation (8).
4. Conclusions
In the previous studies, the data are analyzed separately in three groups according to the dimensions and not treated as a whole. And only the relationship between NPR and is derived for each fixed dimension . As a consequence, there are diverse formulas for different dimensions. In the study, we further find the relationship between NPR and for each fixed . Finally, based on the two relationships and all the data treated as a whole, the relationship among the three is obtained and described in Equation (8), which shows the NPR scales with the (the dimension) power of the product of the square root of the diffusion coefficient (coupling or gap junction ) and a constant which might depend on the strength of the pacing cells. And then, the accuracy of the formula is evaluated by the goodness of fit (), inferring that the formula fits the data well. Moreover, a 1D strand model and a 2D tissue model are developed to derive more data to verify the formula, and a positive result is received.
In conclusion, based on the data from the reference and the simulation, the formula is proposed and validated to describe the relationship among the NPR , the dimension , and the diffusion coefficient (coupling or gap junction ). And only one formula is required for all the data in the three dimensions in the study, while three formulas are needed according to the previous studies. Moreover, especially for the 3D simulation, there are always 106~107 cells in a volume. It will cost a lot of time only for a single simulation period. Therefore, much time is spent on finding the suitable NPR. However, according to the results in the work, formula (8) could be determined in advance based on the 1D and 2D data. And then, the NPR for the 3D simulation could be predicted by the determined formula. As a result, much time is saved.
On the other hand, more simulation and experimental data are expected to verify the formula; what is more, further investigation is necessary to probe the biological meaning of the parameters and .
Data Availability
The data supporting this research are from previously reported studies, which have been cited.
Conflicts of Interest
The authors declare that they have no competing interests.
Authors’ Contributions
All authors conceived the study, carried out the proofs, and approved the final manuscript.
Acknowledgments
The authors are grateful to the referees for providing the comments and suggestions which lead to improvements of the paper. This work is supported in part by the National Natural Science Foundation of China under Grant 61701135, in part by the China Postdoctoral Science Foundation under Grant 2018M630342, in part by the Natural Science Foundation of Heilongjiang Province under Grant QC2018075, and in part by the Fundamental Research Funds for the Central Universities under Grant 3072020CFQ0601 and 3072020CF0604.