The Quantitative Relationship among the Number of the Pacing Cells Required, the Dimension, and the Diffusion Coefficient

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 i, and the diffusion coefficient D (electrical coupling or gap junction G). The relationship between NPR and G has been investigated in different dimensions, respectively. That is, for each fixed i, there is a formula to describe the relationship between NPR and G; and three formulas are required for the three dimensions. However, there is not a universal expression to describe the relationship among NPR, G, and i 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, i, and G; and then, employing the cftool in MATLAB, the explicit formulas are derived for different cases. In addition, the goodness of fit (R2) 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, i, and the G (D) could be described by a universal formula, where the NPR scales with the i (the dimension) power of the product of the square root of G (D) and a constant b which is dependent on the strength of the pacing cells and so on.


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][4][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][9][10][11][12].
As a consequence, the biological pacemaker has been attracting the attention of the researchers to overcome the disadvantages of the electronic devices [13][14][15][16]. One of the popular strategies is to create biological pacemakers in the ventricle [17][18][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][23][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 i, there is a formula to describe the relationship between NPR and the gap junction G, 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 G, and dimension i 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, i, and G. 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 I K1 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.

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.  [31], which is shown in Equation (1).
where I ion = I Na + I K1 + I to + I Kr + I Ks + I CaL + I NaK + I NaCa + I pK + I pCa + I bCa + I bNa : ð2Þ In Equation (1), V is the transmembrane potential; I ion is the sum of all the transmembrane ion currents; I stim is the external stimulus current and is set to 0 in the study; C m is membrane capacitance per unit surface area; I K1 is the inward-rectifier K + current; I Ks , I Kr , and I to are the outward slow, rapid, and transient and rectifier potassium currents, respectively.
In the simulation, I K1 is depressed to 0.05 nS/pF to obtain the robust pacing cells. 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 D is the diffusion coefficient with the normal value 15.4 mm 2 /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.

The MATLAB Tool.
The cftool in MATLAB is adopted for the fitting in the manuscript.

Results and Discussion
In this section, our aim is to derive a formula to fit the relationship among the NPR, the electrical coupling D (gap junction G), and the dimension i: 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, N i ði = 1, 2, 3Þ is approximately a geometric sequence. The details are shown in Tables 4 and 5.   2 BioMed Research International      Tables 4 and 5, N 2 /N 1 is close to N 3 /N 2 , which infers that for a fixed G, N i ði = 1, 2, 3Þ might be geometric sequences corresponding to i, and N i+1 / N i is the common ratio qðGÞ. As a consequence, the expression of N i of i can be described as where mðGÞ and qðGÞ are functions of G, and i is the dimension.
What is more, the average values of qðGÞ for each G in different cases are computed according to the values of N i+1 /N i 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 q 780 / q 125 is close to Therefore, the expression of N i about i can be rewritten as Equation (5).
According to [28][29][30], for a fixed i, N i is proportional to the i power of the square root of the gap junction (coupling conductance) in the i dimensional tissue, and therefore, N i could also be expressed as follows: where f ðiÞ is a function of i, and G is the gap junction (coupling) conductance. Considering Equations (5) and (6) together, for ∀G and ∀i, there is Therefore, mðGÞ must be a constant, set as a. And the formula of Nði, GÞ could be described as Equation (8).
where a and b are the undetermined coefficients.
In the previous studies [28][29][30], only the relationship between NPR and G for each fixed dimension i is found. As   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 Nði, GÞ 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 Nði, GÞ 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 Nði, GÞ, 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.

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 G, the formula is an exponential function of i, and it could be transformed logarithmically to derive a linear function. And then, the goodness of fit (R 2 ), 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 G = 125 nS and G = 780 nS are shown in Figure 6. And the 6 logarithmic original N i in Table 1 for case "EAD/normal" are drawn in pink for G = 125 nS and in red for G = 780 nS, respectively. Figure 6 illustrates that the formula fits the original values well intuitively. Furthermore, R 2 is computed to evaluate the fitting quantitatively.
The closer R 2 is to 1, the better the fitting is. And R 2 values are 0.9964 and 1.0000 for G = 125 nS and for G = 780 nS, 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) The goodness of fit (R 2 ) 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 R 2 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 5 BioMed Research International myocytes by depressing I K1 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 D are shown in Table 9.
Next, the relationship between gap junction G and diffusion coefficient D is introduced first. The diffusion throughout     BioMed Research International a coupled tissue depends on the presence of the gap junction G. And D 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, I stim 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.

BioMed Research International
Comparing Equation (19) with Equation (20), it is found that G gap is corresponding to D/0:15 2 . That is, except for a coefficient, G gap and D are essentially the same. As a consequence, the data about NPR, D, and i 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 (Nði, DÞ), the diffusion coefficient D, and the dimension i, 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 i = 1, the equation could be rewritten as Based on formula (22), the ffiffiffi ffi D p -NðDÞ graph is shown in Figure 11 together with the corresponding original values in Table 9.
And, when i = 2, The graph of function (23) is shown in Figure 12 together with the corresponding original values in Table 9.
The goodness of fit (R 2 ) 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 G in references [28][29][30] and the relationship between NPR   Table 9.  Table 9 for the 1D pacemaker.  Table 9 for the 2D pacemaker. 9 BioMed Research International and i found in the study, together with the data from reference [30] and the simulation in the paper, the relationship among the NPR (Nði, GÞ), the dimension i, and the gap junction G (diffusion coefficient D) is derived and evaluated, which satisfies the two-variable function in Equation (8).

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 G is derived for each fixed dimension i. As a consequence, there are diverse formulas for different dimensions. In the study, we further find the relationship between NPR and i for each fixed G. 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 i (the dimension) power of the product of the square root of the diffusion coefficient D (coupling or gap junction G) and a constant b which might depend on the strength of the pacing cells. And then, the accuracy of the formula is evaluated by the goodness of fit (R 2 ), 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 Nði, GÞ = a * ðb ffiffiffi ffi G p Þ i is proposed and validated to describe the relationship among the NPR Nði, GÞ, the dimension i, and the diffusion coefficient D (coupling or gap junction G). 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 10 6~1 0 7 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 a and b.

Data Availability
The data supporting this research are from previously reported studies, which have been cited.