A Computer Simulation Research of Two Types of Cardiac Physiological Pacing

: This manuscript adopted the cardiac modeling and simulation method to study the problems of physiological pacing in clinical application. A multiscale rabbit ventricular electrophysiological model was constructed. We simulated His-bundle pacing (HBP) treatment for left bundle branch block (LBBB) and atrioventricular block (AVB), and left bundle branch pacing (LBBP) treatment for LBBB by setting various moments of the stimulus. The synthetic ECGs and detailed electrical activities were analyzed. Our electrophysiological model accurately simulated the normal state, HBP, and LBBP. The synthetic ECG showed that QRS duration was narrowed by 30% after HBP correction for LBBB. For LBBB correction with LBBP, the synthetic ECGs of LBBP starting before 30 ms (if the end of atrial excitation is set as 0 ms) presented right bundle branch block (RBBB), and those of LBBP starting at 30–38 ms were synchronous, while those of LBBP starting after 42 ms possessed LBBB morphologies. The best pacing results were obtained when LBBP started at 34 ms. This manuscript veriﬁed the feasibility of the constructed ventricular model, and studied the physiological pacing mechanism. The results showed that HBP realized correction for AVB and high LBBB. The performance of LBBP can be improved by applying the stimulus within a speciﬁc period of time (0–8 ms) after atrial excitation.


Introduction
Cardiac resynchronization therapy (CRT) is an important clinical treatment for many bradyarrhythmias and heart failure. Several forms of delivering CRT have been applied in recent years. The emergence of biventricular pacing technology introduced treatment plans with pacemaker in the form of CRT. However, it is unresponsive to approximately 40% of treated patients [1,2]. Moreover, it does not constitute left and right ventricular synchronization in the true physiological sense. Physiological pacing is a new effective method in CRT, which refers to pacing directly in the conduction system. In 2000, Deshmukh first applied His-bundle pacing (HBP) for the treatment of patients with atrial fibrillation in combination with atrial ventricular nodal ablation [3]. HBP was gradually applied in the clinic and its clinical validity and feasibility have been demonstrated in many studies [4]. However, the difficulty of accurate positioning and high pacing thresholds of HBP limited its widespread application in clinical practice. In 2017, left bundle branch pacing (LBBP) was first proposed and practiced [5]. LBBP uses a relatively simple implantation technology and more stable and reliable pacing parameters, and constitutes an alternative physiological pacing method [6,7].
Indications for physiological pacing include the atrioventricular block (AVB) and bundle branch block. Although LBBP has become a new popular treatment, it still has some limitations compared with HBP. For example, the conduction through the right bundle was abandoned, the excitation of the left ventricle is often too much earlier than that of the right ventricle. Additionally, the right bundle branch block (RBBB) often occurs following the use of LBBP for the correction of both the left bundle branch block (LBBB) [5,8] and AVB [9], which may result in right heart insufficiency in the long run.
Therefore, the key problems that need to be solved in the clinical application of physiological pacing include (a) the detailed internal mechanism of HBP and LBBP, (b) the setting of the pacing parameters to achieve biventricular synchronization in LBBP, such as release time of stimulus pulse. In clinical practice, the performance of physiological pacing can only be evaluated by clinical features, such as the electrocardiogram (ECG), ejection fraction, and others. The internal principles of electrical activities can hardly be explored and clarified in detail through clinical experiments. These years, computer modeling has become an effective way to study cardiac electrophysiology and has been applied in many studies [10,11]. This method can simulate the process of action potential conduction on a heart model that accounts for the electrophysiological characteristics of cells and tissues. Recently, this method has begun to be used in the study of physiological pacing [12].
This manuscript adopted the method of cardiac modeling and simulation to conduct research related to pathological problems in the clinical application of physiological pacing. First, we constructed a detailed, multiscale, rabbit ventricular electrophysiological model, which contained myocardial segmentation, cellular heterogeneity, and His-Purkinje system. We compared the synthetic ECG with the real ECG dataset to verify the correctness of the model. The HBP treatment for LBBB and AVB was then simulated, and the internal mechanism was explored, which also further verified the feasibility of the modeling simulation method for physiological pacing research. Finally, LBBP treatment for LBBB was simulated, we investigated the different effects of the release time of stimulus pulse, and determined how to improve pacing efficiency of LBBP.

Rabbit Ventricular Electrophysiological Model with His-Purkinje System
The electrophysiological model used in this manuscript is an optimized rabbit ventricular model. We completed the optimization based on an anatomically detailed finiteelement rabbit ventricular model that was created from magnetic resonance images. The model consisted of 82,619 nodes and 431,990 tetrahedral elements [13]. The myocardium was separated into three layers with equal thicknesses: endocardium, midmyocardium, and epicardium, according to the transmural distance. The Mahajan cell models [14] were added on three layers of myocardial tissue, with most of the parameters set on default, while several parameters were adjusted to ensure the correct sequence during repolarization, which can be reflected in the action potential (AP) morphology of each myocardial layer and the polarity of the T wave. We changed the slow inactivating transient outward K+ current (I to,s ) conductance and slow delayed rectifier K+ current (I Ks ) conductance according to the ratio of those in [15]: g to,s = 0.04 mS/µF for epicardial and M cells, g to,s = 0.01 mS/µF for endocardial cells; and g Ks = 0.1386 mS/µF for endocardial and epicardial cells, and g Ks = 0.0347 mS/µF for M cells. We simulated the electrical activation with the bidomain model of this rabbit ventricle that is commonly used in cardiac electrophysiological simulation presently. The bidomain model considers tissue as a two-phase medium along with extracellular and intracellular spaces [16]. We introduced anisotropic conduction based on the microstructural fiber orientations contained in the model. Extracellular and intracellular conductivities were set according to the study of Clerc L [17]. Specifically, a conductivity of 0.62 S/m was set along the fibers and 0.24 S/m perpendicular to the fibers in the extracellular space, while the corresponding conductivities in the intracellular space were 0.174 S/m and 0.019 S/m, respectively. We then multiplied all these ratio values by 0.3 times to make the conduction velocity of the simulated myocardium consistent with that of the real rabbit ventricular myocardium.
A detailed His-Purkinje system was integrated in the ventricular model. The His-Purkinje system model was generated according to the method proposed by Costabal [18], which had been applied in many modeling studies [19,20]. In this method, the twodimensional fractal tree network was first generated, and was then projected on the endocardium by a two-step projection method to form the Purkinje network on the ventricular model. To conduct the physiological pacing simulations more efficiently, the following details were noted in the construction of the His-Purkinje system in this manuscript: We added an extra branch on the top of the Purkinje network to represent the His-bundle. Previous studies had indicated that the left and right bundle branches have separated within the His-bundle [21]. The longitudinal dissociation is the reason HBP can correct the bundle branch block that occurs within the His-bundle [22]. Therefore, we constructed the structure in which the left and right bundle branches were dissociated within the His-bundle. There are left and right bundle branches, and the left bundle branch divides into the left anterior and the left posterior fascicles. The His-Purkinje system model we generated contained 369 branches and 146 terminals, and its morphology was consistent with its true physiological structure. Figure 1 illustrates the optimized rabbit ventricular model with the His-Purkinje system model. The gray part and the white network represent the myocardium and the His-Purkinje system model, respectively. The endocardium was connected with the Purkinje network model by Purkinje-muscle junctions (PMJ). Based on this ventricular electrophysiological model, we simulated the cardiac physiological rhythm to verify the feasibility of the model and the proposed method. Given that the model did not include the atrium, to simulate atrioventricular sequential pacing, it was hypothesized that the His-bundle excitation occurred after the atrial and atrioventricular node. According to ECG data from rabbits, the duration of the PR segment (the time between the onset of the P wave and the onset of the QRS complex) is approximately 40 ms. Therefore, in this study, the end of the P wave was set as 0 ms, and the time of slow conduction in the hypothesized atrioventricular node lasted from 0 to 30 ms. At 30 ms, the proximal end of His-bundle initiated the excitation. Thus, the His-Purkinje system was excited from 30 to 40 ms, and the ventricular myocardium would begin its excitation at approximately 40 ms, which corresponded to the onset of the Q wave.

Synthetic ECG
Synthetic ECG data were generated for these simulations with the Cardiac Arrhythmia Research Package (CARP) [23]. When the bidomain was formulated, extracellular current flow and potentials were computed everywhere in the entire domain. This approach was the most accurate and was considered the gold standard for ECG modeling [24]. We got limb leads and chest leads recordings by setting electrodes on the heart model according to the lead-axis of 12 leads.

Simulation of HBP and LBBP
We conducted all the simulations within CARP. The time step of the simulation was set to 0.005 ms. The bundle block was implemented based on the adjustment of the conductivity of a specific branch to 1/2000 of that in physiological conditions. This assured that the conduction would not be activated through the branch. Two LBBB types were simulated in the ventricular model. In LBBB-1, the block was set in the left bundle branch (LBB) part within the His-bundle. When the excitation was conducted through His-bundle, LBBB appeared, as shown in Figure 2a. In LBBB-2, the block was set in the left main bundle branch, and resulted in a typical LBBB. AVB was also simulated in the ventricular model by setting the block in the distal His-bundle so that atrial excitation could not conduct to the ventricle, as shown in Figure 2b. The simulated physiological excitation initiated from the proximal end of His-bundle started at 30 ms in all simulations. The method of pacing simulation involved the application of a stimulus that lasted 0.5 ms that was strong enough to stimulate the branch and myocardium. The pacing simulation was non-selective HBP (NSHBP) or non-selective LBBP. This meant that both the bundle branch and local myocardium were stimulated.

HBP Simulations
We performed HBP to correct LBBB-1 and AVB, respectively. This was equivalent to the restoration of the function of His-bundle. To conduct HBP, we applied a stimulus that started at 31 ms at the distal end of the His-bundle, which was beyond the bundle block position, as shown in Figure 2c,d. ECG simulations were carried out for LBBB, LBBB correction with HBP, and AVB correction with HBP.

LBBP Simulations
In this part, we simulated LBBP to explore the effects of pacing moments on ventricular synchronicity under LBBB conditions. LBBP was implemented beyond the block area of LBBB-2 by applying a stimulus at the distal end of the left main bundle branch, where the left anterior and left posterior fascicles began to separate, as shown in Figure 2e. Given that the His-bundle and the right bundle branch (RBB) were still functionally normal, the physiological excitation initiated from the proximal His-bundle could be conducted to the right ventricle. Therefore, the relative time between the LBBP stimulus and the excitation of the proximal His-bundle could have a great influence on the performance of the LBBP. We simulated nine different situations in which the LBBP stimulus started at 18 ms, 22 ms, 26 ms, 30 ms, 34 ms, 38 ms, 42 ms, 46 ms, and at 50 ms, respectively, while the proximal His-bundle excitations all started at 30 ms. Therefore, the relative time between the LBBP stimulus and the excitation of the proximal His-bundle varied from −12 ms to +20 ms. The conduction process and ECG in each situation were compared and analyzed to explore the best pacing moment for the use of LBBP for the correction of the LBBB. Figure 3a shows the simulated APs of the cells from the endocardium, epicardium, and midmyocardium, recorded under steady state. The action potential duration (APD) of the midmyocardial cells was the longest, and the APD of the endocardial cells was slightly longer than that of epicardial cells (208 ms in the midmyocardial cells vs. 196 ms in the endocardial cells, and 187 ms in the epicardial cells).

APs and ECG of Ventricular Model
The synthetic ECG results in normal ventricles are shown in Figure 3b. Given that the leads I, II, V1, V2, V5, and V6 have obvious characteristics, we chose these leads for analyses. The duration of the QRS complex was approximately 38 ms, and was consistent with that of a real rabbit ECG. The morphology of the synthetic ECG was also similar to that of the real ECG. This indicates that the ventricular electrophysiological model constructed in this manuscript can reflect the normal ventricular electrophysiological activity.

HBP
The clinical ECG data before and after treatment of LBBB with HBP are both shown in Figure 4a,b for the comparison with the simulated results [25]. The reference study implemented NSHBP on 12 patients with LBBB and corrected their LBBB with a narrowing of QRS duration by 30%. The clinical ECG data were achieved from a typical patient with ischemic cardiomyopathy.  In the simulation of HBP correcting LBBB-1, the His-bundle and the basal segment of the ventricular septum were both stimulated simultaneously. The LBB part within Hisbundle remained blocked, and the activation of the RBB was caused only by the excitation from His-bundle. The synthetic ECG of LBBB-1 is shown in Figure 4c. The duration of the QRS complex was 55 ms in lead II, and the left ventricular activation time (LVAT) was 32 ms. The R waves in leads I, II, V5, and V6, were wide and large. Leads V1 and V2 were rS-shaped, with a much smaller R-wave and a deeper and wider S-wave compared with those in normal states. The T waves were inverted in leads I, II, V5, and V6, and were upright in leads V1 and V2. Figure 4d shows the conduction process of the electrical activity in LBBB at several specific moments that corresponded to the synthetic ECG characteristics. We used the QRS complex in lead II as a reference for detailed explanations. The activation sequences at 48.8 ms illustrated that the excitation of the PMJ was corresponding to the onset of the R-wave. The activation sequences at 62 ms showed that the wavefront in the right ventricle had fused and began to conduct to the left ventricle, which corresponded to the first R-wave peak. At 68.6 ms, the wavefront had excited the Purkinje network in the left ventricle and the PMJ. This time point corresponded to the onset of the second R-wave. Finally, at 77 ms (R-wave peak), the wavefront in the left ventricle fused.
After the HBP was implemented, the LBB restored its conduction function, and the left ventricle was excited normally. Figure 4e shows the QRS complex of the synthetic ECG under HBP (blue line) compared with the synthetic ECG of the normal state from Figure 3 (red line). The duration of the QRS complex was 40 ms, narrowed by 28%, similar to that of the normal state. In the real ECG of the normal state, the Q wave is a reflection of the earliest excitation that occurs on the left side of the septum. Most of the myocardium depolarize synchronously, and correspond to the R-wave. Part of the left ventricular basal myocardium eventually depolarize, and thus generate the S-wave [26]. The simulated activation sequences of HBP are shown in Figure 4f. There was no Q wave in the synthetic ECG because the basal segment of the septum was excited first as soon as the NSHBP began. The pacing stimulus and the onset of the QRS complex were not discrete on the synthetic ECG because the excitation was conducted both through the His-Purkinje system and the myocardium. At 48.8 ms, the PMJ in the two ventricles were excited that indicated the LBB restored conduction function. At 58 ms, the two ventricles were excited synchronously. This corresponded to the R-wave peak in the synthetic ECG of lead II and lead V5, meaning the LVAT was 28 ms. Given that the pacing excited the basal segment of ventricular septum at the earliest, the left ventricular basal myocardium and most parts of the ventricles depolarized at the same time at 62 ms, and the amplitude of the corresponding S-wave peak was not as large as that in the normal state.
The synthetic ECG of the HBP correcting LBBB-1 corresponded with the real ECG both statistically and morphologically. The proportion of QRS reduction was similar (28% vs. 30%). The R waves had physiological shapes. In the LBBB results, the polarities of the T waves were opposite to those in normal states. After the completion of HBP, the T wave polarities were corrected. The results indicated that the simulation method reproduced the bundle branch block and conduction system pacing well. This further proved the validity of our model and method.
The synthetic ECG of HBP correcting AVB is shown in Figure 4g, which is also very similar to that in the normal state. HBP excited the conduction system at the distal Hisbundle, and resynchronized the ventricles. The simulated activation sequences were nearly the same as those in HBP used for the correction of LBBB given that the pacing positions of these two situations were very close and the conduction paths of the wavefront were basically the same.

LBBB Correction with LBBP
In this case, although LBBB-2 occurred in the left main bundle branch, it still yielded a similar synthetic ECG as that in Figure 4c. With the existence of a normal beat from His-bundle, we performed LBBP at different moments. The synthetic ECGs under LBBP varied morphologically. Given that the variations on the QRS complex were more obvious than the variations on the T wave, we used the QRS complex for analysis. Figure 5 indicates respective QRS complexes and the QRS duration for each specific moment, compared with the referenced synthetic ECG of the normal state. The duration of the QRS complexes of these simulations were 21%, 16%, 13%, 5%, 5%, 8%, 47%, 58%, 63% longer than that of the synthetic ECG of the normal state, respectively. The synthetic ECGs of LBBP starting at 18 ms, 22 ms, and 26 ms had RBBB morphologies, and those of LBBP starting at 30 ms, 34 ms, and 38 ms appeared to be basically synchronous, while those of LBBP starting at 42 ms, 46 ms, and 50 ms possessed LBBB morphologies. The best pacing results were obtained when LBBP started at 34 ms. Morphologically, the synthetic ECG of this situation was extremely consistent with the normal synthetic ECG. We present the simulated activation sequences of three typical situations (LBBP starting at 18 ms, 34 ms, 50 ms), as shown in Figure 6. Figure 6a shows the results of the situation in which the LBBP stimulus started at 18 ms, i.e., at 12 ms earlier than the proximal His-bundle excitation. The right ventricle was excited after the left ventricle and resulted in an RBBB. Regarding the best pacing situation, the pacing was 4 ms later than the proximal His-bundle excitation, as shown in Figure 6b. The entire conduction process of the wavefront was very similar to that in the normal state, which indicated that LBBP corrected LBBB successfully and led to a physiological-like rhythm. Figure 6c shows the results of the situation in which the LBBP stimulus started at 50 ms, i.e., at 20 ms later than the proximal His-bundle excitation. The left ventricle was excited after the right ventricle. This caused an LBBB. Therefore, LBBP failed to correct the LBBB in this situation.

Model Construction
We constructed a detailed ventricular model that included the His-Purkinje system to conduct the physiological pacing simulation study. In the ventricular conduction system model, the RBB was one branch, and the LBB was divided into the left anterior and the left posterior fascicles. We set up the conduction block by reducing the conductivity of the branch to 0.05% so as to make it lose its normal conduction function. Compared with the model constructed in other studies of physiological pacing [12], our model took into account the longitudinal association in the His-bundle and the electrophysiological heterogeneity of different layers of cardiomyocytes, which had been demonstrated as important in other studies [27]; and we simulated the electrical activation with the bidomain model which was more accurate. The relationship between the APs of the three layers of myocardial cells was consistent with that of previous studies [15]. The synthetic ECG obtained by the model in this manuscript was very close to the real ECG, which indicated that the process of cardiac electrical activities shown by the model matched the actual situation.

HBP Simulations
As it can be observed in Figure 4, the synthetic ECGs before and after the HBP in the case of LBBB were basically consistent with the clinically measured results, which validated the model. The QRS duration and LVAT were shortened after the HBP, and the morphology of LBBB was eliminated. On this basis, with the advantage of cardiac modeling and simulation, this manuscript further analyzed the cardiac electrical activity process before and after HBP, and revealed the mechanism of HBP responsible for the correction of LBBB. As it can be observed from the activation sequences, the two ventricles were resynchronized. The LBBB was corrected with HBP because the block was within the His-bundle. In the case of AVB, the synthetic ECG of HBP was very similar to that in the normal state as well. Our pacing method was NSHBP, which meant that the stimulus was conducted through both the His-Purkinje system and the myocardium. Therefore, in the corresponding synthetic ECG, the pacing stimulus and the onset of the QRS complex were not discrete. This is consistent with the conclusions of the clinical HBP study [8].

LBBP Simulations
Compared with HBP, the results of the LBBP were more complicated. The pacing moment had an obvious impact on the performance of the stimulated LBBP. When the pacing occurred much earlier or later than the excitation of His-bundle, RBBB or LBBB would be induced, respectively. Only when the pacing was applied later than the excitation of His-bundle within a certain period of time, the LBBP corrected LBBB successfully and the ventricle achieved its resynchronization. The reasons behind this result are as follows: in the situation in which the best LBBP results were obtained, the pacing stimulus was applied 4 ms later than the His-bundle excitation at the distal end of LBB. In fact, in the simulation of the normal state, the His-bundle excitation required exactly 4 ms to achieve conduction to the distal end of LBB. LBBB treated by LBBP required this specific time to fuse the activation wavefront with normal conduction in the His-bundle and RBB to form a normal QRS morphology. The relative time between the atrial excitation and the pacing can be considered to be the atrioventricular delay (AVD) of the pacemaker. Therefore, a proper AVD can help LBBP avoid the occurrence of the RBBB and attain good performance. In some clinical LBBP studies, researchers had mentioned that the adjustment of AVD could eliminate the RBBB morphology in the ECG results [28]. Consistent with our findings, Strocchi M et al. proved that AVD can be optimized to improve right ventricular activation times [12]. Yet we explained the reason behind the result, explored the more detailed conduction mechanism of LBBP, and achieved the conclusion with a small-scale animal heart model. Thus, our method was more high-efficiency. Therefore, this cardiac modeling and simulation can be used to optimize the study of the stimulus release time, and can show the detailed excitation conduction process in various situations to analyze and reveal the mechanism of the influences of the pacing time on the pacing effects.

Limitations
There are some limitations associated with our study. Considering the limitations of computing resources, this manuscript used a rabbit ventricular model to improve the simulation efficiency. Also, the model relied on some assumptions which were not completely physiologic. To address these issues, in a subsequent study, we will construct a human heart model that includes a personalized His-Purkinje system, and embed it in a torso model to include the chest leads, and incorporate the atrial model. We will also construct various types of the cardiac conduction system in accordance with anatomical information to provide the theoretical basis for clinical-relevant LBBP research.

Conclusions
Based on the construction of an elaborate rabbit whole ventricular model, this manuscript proposed a modeling and simulation method to study the pacing mechanism and parameter optimization methods of HBP and LBBP, and verified its feasibility. The simulation results showed that non-selective HBP can be used to treat AVB and high LBBB (LBBB-1). The therapeutic effect of LBBP was related to the pacing time. To obtain the best performance when applying LBBP to treat LBBB, the stimulus needs to be delivered within a specific period of time (0-8 ms) after the excitation of His-bundle.  Data Availability Statement: Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.