A Two-Stage Purkinje Network for More Accurate ECG Representations

In this manuscript we propose a method to generate Purkinje networks that are anatomically and physiologically plausible, for use with in-silico modeling. Purkinje networks play a fundamental role in shaping cardiac electrical activation patterns and their corresponding clinical electrocardiograms (ECGs). Despite a known variability in ventricular activation sequences, certain sites of early activation within the left and right ventricles have been identied in the literature for normal electrical excitation patterns. Nevertheless, in-vivo imaging of Purkinje networks cannot at present yield detailed information on their structure, so there is a genuine need for in-silico models that can construct Purkinje networks that are both anatomically and physiologically plausible, in particular networks that can exhibit correctly situated early activation sites in the ventricles. Of special interest to this manuscript is the method of representation of Purkinje networks by line-like electrical elements that are generated by means of a fractal-tree algorithm (1–3) to overlay the irregular endocardial surfaces. A known drawback to a direct implementation of this approach in complex geometries relates to its incorrect modeling of clinically observed ECGs and electrical activation sequences for a human heart (4). Our aim was thus to correct this deciency by generating Purkinje networks that leverage a pre-knowledge of the location of early activation sites. At every such site we rst generate a Purkinje sub-network. These sub-networks are linked together and to the bundle of His, setting up our rst stage of the Purkinje network. Subsequently, we spawn a second stage to the Purkinje network from one or more tips of any given sub-network, to cover the full endocardial surface with Purkinje elements. Our resulting activation sequences and ECGs compare favorably to those of a population of 39 healthy male individuals (the PTB diagnostic database), and our corresponding mechanical markers of cardiac function also match well with the literature.


Introduction
Despite extraordinary modern day advances in drug development, surgical strategies and cardiac therapies, cardiovascular disease (CVD) accounts for 31% of global mortality (5), with the American Medical Association (AMA) reporting almost 18 million CVD-related deaths annually (6). According to the AMA, heart disease remains the number one cause of death in the United States, and accounts for the country's highest direct health expenditure. Moreover, the largest number of surgical procedures in the United States are on the cardiovascular system (6).
As a result, research on CVDs has been moving forward apace, and towards in-silico modeling and simulation of cardiac electromechanics for the past couple of decades. The cardiovascular system is a multiscale system that operates from the molecular to the cellular to the tissue level to de ning the whole organ (7). Monitoring changes in the smaller levels could reveal how they affect the heart's structure and function and give us the ability to better quantify how heart disease impacts its functionality (8). Thus, with the aim to improve cardiac healthcare, multiscale modeling could be pursued to gain insight into the causes and consequences of CVDs (7,9,10), monitor drug-induced effects on cardiac electromechanics (3,11), optimize treatment options and simulate surgeries to nd the best course of operation (12)(13)(14).
However, these models are only as valuable as their ability to reproduce clinical data accurately, and to correctly predict clinically relevant outcomes.
Computational multiscale models have initially mostly focused on modeling only a subsystem of interest (usually the left ventricle) (9,15,16) but have more recently evolved to modeling both ventricles (17)(18)(19) and then to modeling all 4 chambers (4,20,21). Currently, more models also include valves, papillary muscles and major vessels (aortic arch, pulmonary arteries and superior vena cava) (7). Generally, the reservation modelers have to create and use whole heart models relates to the high computational cost of the resulting simulations, as well as the complexity in image segmentation, discretization and meshing (4). For instance, many computational models have excluded the Purkinje network from their geometry altogether, choosing to apply excitation stimuli to ventricular tissue instead (11,22), since in-vivo imaging of the networks cannot at present yield such detailed information on structure (23). This overlooking of the Purkinje network could cause errors in the model's performance, since Purkinje ber cells are fastconducting cells that transport electricity at about 6 times the rate of ventricular muscle bers (24).
Purkinje bers thus play a fundamental role in synchronizing the contraction of the left and right ventricles and in shaping cardiac electrical activation patterns and their corresponding clinical electrocardiograms (ECGs) (23). Moreover, physical defects to the Purkinje bers lie at the root of various cardiac diseases, e.g. Bundle Branch Blocks (24), so that their inclusion in computational models is generally important to capture diseased cardiac behavior in relation to its underlying structural causes.
Since their invention, ECGs are still the most widely used non-invasive method to assess cardiac function and to investigate cardiac disorders (25). Despite the known natural variability to ventricular activation sequences across a healthy population, which translates to a variability observed in the QRS complexes of human ECGs, certain sites of early activation within the left and right ventricles have been identi ed in the literature, which de ne normal electrical excitation patterns (25). Thus, when interested in building an accurate multiscale cardiac model that accounts for electrophysiology within anatomically accurate cardiac domains, it becomes necessary to ensure that the resulting ECGs from the model match those clinical, both in health and in disease.
One such model is the state-of-the-art multiscale Living Human Heart Model (LHHM), recently developed by Simulia as part of a translational research initiative that it launched and now includes over 90 organizations across the spectrum of research and educational institutes, industrial partners, clinicians and regulatory bodies (4). The LHHM is an electrically excitable and deformable, anatomically accurate, 4-chamber heart model based on the magnetic resonance imaging of a healthy middle-aged Caucasian male (26). Presently, there remains in the LHHM's electrophysiological model a drawback that relates to its inability to correctly model clinically observed electrical activation sequences and their resulting ECGs for a normal heart. The aim of this paper is thus to model on the LHHM ventricular activation via a new two-stage Purkinje network model, whose ECG predictions are validated against clinically derived human activation data. This is achieved through a generation of a two-stage Purkinje network that depolarizes certain pre-identi ed areas of early activation on the ventricular endocardium ahead of the remaining ventricular tissue. This two-stage Purkinje network approach is integrated into the LHHM electromechanical solver, and pseudo-ECGs produced from this modi ed model are compared to a sample of healthy clinical ECGs and to the pseudo-ECGs produced by the out-of-the-box LHHM. Electrical and mechanical comparisons are also made for a population of healthy individuals.

Two-stage Purkinje network
We herein focus on the Purkinje network model proposed by Costabal et al. (23), which is generated by a fractal tree algorithm and was used to produce the LHHM Purkinje network. Its python code is freely available at https://github.com/fsahli/fractal-tree. The algorithm works by taking as input the model endocardium object, the coordinates of the Purkinje node of origin on the ventricular endocardium and the direction of propagation. Unique line elements are then automatically generated in an iterative tree branch-like pattern that are projected on the endocardium surface. The length, angle and curvature of the line elements are controllable. The complexity and superiority of this algorithm lies in its ability to project the network nodes on the irregular endocardium surface, including the protruding papillary muscles. The irregularities of the ventricular surface and the fact that the papillary muscles are known to be covered in Purkinje bers have usually been overlooked by typical cardiac models unable to capture such complexity.
Even though there exists a notable level of variability from person to person in electrical activation sequences, experimental studies have highlighted certain similarities across human datasets. These studies include ex-vivo microelectrode recordings (27), in-vivo electromechanical wave imaging (28) and non-invasive electrocardiographic imaging (29). Based on these studies, Cardone-Noott, et al. (25) identi ed 7 root points in the ventricles from which excitation propagates; 4 in the left ventricle and 3 in the right ventricle. They tested the effect of the activation map on their heart-torso model ECGs by assuming the sites of earliest conduction to be the Purkinje-muscle junction sites and stimulating the endocardium there. Further con gurations were done by varying passive tissue conductivities and tissue activation speed.
In this work we set off to take advantage of the LHHM and generate a Purkinje network with the primary aim of reproducing clinical ECGs of a healthy heart. The endocardium structure of the LHHM ventricles was extracted for use in the Purkinje generation. For each of the ventricles, we modi ed the Costabal et al. (23) fractal tree algorithm to generate a network that primarily ties to the endocardium at the root points of early activation before spreading to the remaining ventricular tissue. This was done for each ventricle by generating Purkinje ber networks of limited size on the endocardium at these early activation sites indicated on Figure 1(a). We will refer to these disconnected Purkinje branches as stage 1.
Stage 2 of the network was generated from the end of a single stage 1 branch in each ventricle to cover the remaining endocardium. It is possible to generate stage 2 from the end of each (rather than a single) stage 1 site but was deemed unnecessary to reach our objective of reproducing physiologically plausible ECGs. Through several trials and adjustments, the fractal tree algorithm parameters we used were a branch length of 6 mm, a branch angle of 0.3 rad, and a repulsion parameter of 0.1. The generated networks in the left and right ventricles shown in Figure 1(b) are then linked through the bundle of HIS which is directly connected to each branch in stage 1. The excitation of the atrioventricular (AV) node will thus travel through the bundle of HIS to activate the ventricular endocardium under the stage 1 branches, before propagating onward to the remaining tissue.
Integration into LHHM electrophysiology model A major drawback to the pseudo-ECGs produced from the out-of-the-box (control) LHHM was an inverted T-wave that manifested in all six precordial leads. Since differences in repolarization across different regions of the ventricles are what contribute to the shape of the T-wave (30), changes in the repolarization rate of the Purkinje bers translate to changes in the T-wave amplitude and direction. By reducing the refractoriness parameter, and in turn the repolarization, of the Purkinje bers, the T-wave direction was corrected on the pseudo-ECGs of the precordial leads. This reduction is physiologically motivated by a closer inspection of the action potentials of the Purkinje bers and the ventricular muscle in the control LHHM. By default, the ventricular muscle had been assigned a lower rate of repolarization than that of the Purkinje, when in fact the Purkinje bers should have the lower rate of repolarization (31). A reduction we implemented to the refractoriness of the Purkinje bers was thus only to correct their relative action potential durations (32), which resulted in the organ-level consequence of a corrected T-wave.
After adjusting the action potentials, our connected two-stage Purkinje network was then imported and integrated into the LHHM, as a replacement of the original Purkinje network. For the coupling between the Purkinje network and the ventricular myocardium to take effect, a node-to-node tie-condition (33) that connects the Purkinje line elements to the endocardial surfaces was implemented. The electric cardiac cycle was lastly adjusted from 500ms to 600ms to produce a pseudo-ECG more comparable to clinical ECGs than those of the control. The LHHM electric simulation was then run using the all other parameters and conditions as those of the control.

Pseudo-ECG generation
The out-of-the-box LHHM does not have a built-in method to generate pseudo-ECGs, which are approximations to the signal propagating off a patient's body surface as a result of changing cardiac electrical potential. An actual ECG signal would be recorded via electrodes placed at speci ed distances from the heart, and would differ from one person to another. For an in-silico model of the heart that does not include the surrounding body, only approximate ECGs can be constructed, viz. a pseudo-ECG. Figure 2 marks approximate positions for six precordial leads, which we use to record the pseudo-ECG signals of our model, i.e. V1, V2, V3, V4, V5 and V6 (34), which are located where actual precordial leads would be anticipated.
We then utilize an open source software library, Chaste (Cancer, Heart and Soft Tissue Environment), which contains a broad set of numerical solvers and algorithms focused on biological modeling.(35) One such solver is its pseudo-ECG calculator which, by mean of Equation 1, calculates the extracellular potential Φ e anticipated at electrode position ( x', y', z').
In Equation 1, is the diffusion coe cient, is the transmembrane potential and is the distance between the electrode and a point within the heart tissue, and integration is performed over the myocardial volume. (36) Hence, given the electrophysiological results produced by the LHHM solver, activation time maps can be generated and converted to an HDF5 le that are read by Chaste's pseudo-ECG calculator. Moreover, the full mesh of the model is provided to Chaste in the form of an element le, a node le, and a face le, to complete the pseudo-ECG calculation.

Pseudo-ECGs vs. Clinical ECGs
The Physikalisch-Technische Bundesanstalt (PTB) diagnostic database (37) is a database made available by the National Metrology Institute of Germany, which contains a collection of digitized ECGs from healthy and diseased individuals, to be used for research and teaching purposes. From it, we obtained the precordial lead ECGs for 39 healthy male individuals, with ages ranging from 17 to 69 years. Figure 3 plots these 39 normalized precordial lead ECG signals (light grey), against our calculated pseudo-ECGs that have been generated by our two-stage Purkinje network (black). Our proposed model seems to t well with the clinical ECGs for all 6 precordial leads, in terms of trends and timings. The upwards R peak grows larger while the downwards S peak grows smaller from signals V1 to V6, as seen in literature (38).
The lead V1 was further examined for the measurement of the PR, QRS and QT intervals as shown in Figure 4(a). The interval timings are recorded in Table 1 against the clinical averages from literature (38). The PR and QRS intervals are seen to lie within the clinical range at 0.14 and 0.07 seconds respectively. The un-adjusted QT interval is below the clinical range at 0.27 seconds; however, since a normal QT interval depends on the heart rate, calculating a corrected QT interval (QT c ) using Bassett's formula gives a corrected interval of 0.34 seconds, which coincides with the clinical average of under 0.4 seconds. Bassett's formula relates the QT interval to the heart rate according to equation 2 as follows (38), where the RR interval is the time interval between the R peaks in two consecutive cardiac cycles and it is used to calculate the heart rate.
Clinical ranges of ECG amplitudes ( Figure 3) exhibit a very large variation from one person to the next, depending on the size of their ventricular chamber and the exact position of the lead on their torso and how far it is from their heart. (39,40) We therefore only consider the average amplitudes of the normalized QRS complex and of the T-wave, e.g. for lead V4 as de ned in Figure 4(b), for the 39 clinical ECGs in our database. We nd that our modi ed LHHM model (with the two-stage Purkinje network) predicts comparable values to the average of the clinical measurements, as summarized in Table 2. Furthermore, the ratio of the amplitude of normalized QRS complex to that the normalized T-wave also compares well to the clinical average, indicating correct their relative prominence in the ECG. Figure 5 summarizes the six precordial lead ECGs for a single healthy individual chosen at random from the database, plotted against the pseudo-ECGs of the control LHHM and our modi ed LHHM. This gure clari es further how the two-stage Purkinje model produces ECGs of closer similarity to the clinical ECGs than those of the control LHHM. For instance, in V1, both the clinical and modi ed LHHM ECGs show the rSr' morphology while the control LHHM shows a QR morphology (24,39). Furthermore, in the remaining ve leads, the downwards S peak is much more prominent in the clinical and modi ed LHHM curves than in those of the control LHHM.
Activation Time Maps Figure 6 shows an endocardial view of the ventricular activation time maps for (a) the LHHM control and (b) the two-stage Purkinje models. The map indicates activation time at each point in the ventricular mesh. The activation map of the LHHM control shows, problematically, that electrical stimulus propagation mostly starts at the apex and spreads upwards towards the base. On the other hand, our two-stage network's consideration of early activation sites indicates that the electrical stimulus actually starts at several different locations across the left and right ventricular walls and septum, then spreads from there. These time maps are generally consistent with clinical measurements reported in the literature (25,(27)(28)(29).
Nevertheless, on closer inspection of the activation maps, we nd that the stage-one Purkinje subnetworks labeled 2 and 6 in Figure 1(a) did not activate their surrounding tissue as early as the other ve did. That is, the myocardial tissue linked to them did not respond to their early activation as quickly as the remaining ve sub-network myocardial neighborhoods, in spite of receiving the early electrical signal. This delay at the two sites appears in the presence of the stage-two Purkinje network. Without stage-two Purkinjes, all seven sub-networks activate their surrounding tissue at roughly the same time. The reason for this delay in tissue activation at sites 2 and 6 relates to a pre-activation of their surrounding tissue from secondary activation sources, which renders them less responsive. Secondary sources include electricity propagated through myocardial tissue from other pre-activated stage-one sites, and through retrograde propagation into stage-two Purkinje bers which then connect to sites 2 and 6. It will be an interesting investigation to consider manners by which secondary sources could be affected by Purkinje network density and its other electrophysiological parameters. For the purposes of this manuscript, however, we note that obtaining clinically comparable pseudo-ECGs could be generated by only adopting ve early activation sites. Indeed, when we modeled a second two-stage Purkinje network that excludes sites 2 and 6, we found essentially the same improvements to the pseudo-ECGs as our network with all seven sites. Conversely, by only locating Purkinje sub-networks at the activation sites (stage-one), without then covering the endocardium by stage-two Purkinje elements, does not yield correct ECGs, though activation maps do appear more consistent with those of the seven sites. As such, we can conclude that while early activation sites have a prominent role in correcting ECGs, they do not need to be fully replicated to reproduce them, and that to obtain correct ECGs a complete coverage of the endocardium by Purkinje bers is advised.
One area where the two-stage Purkinje network can be applied is in patient-speci c, pre-operative Cardiac Resynchronization Therapy (CRT). Despite the successful development and advancement of CRT technology, and its increasing use as a heart failure treatment that synchronizes ventricular activation (41), one third of its cases prove clinically non-responsive. The lack of responsiveness has been linked to multiple factors including the variability of cardiac electrical propagation between patients, and the di culty in determining an optimal placement of the electrodes on the heart (42). Including our two-stage Purkinje network in computational models can help account for patient variability, as it can be easily modi ed to match a patient's ventricular activation pattern. Developing such a patient-speci c model could thus be used to simulate and predict proper CRT lead placement on the ventricles before surgery, and to even determine the likelihood of patient responsiveness to the treatment.

Mechanical Results
To ensure that our modi cations to the electrical LHHM model do not negatively affect ventricular performance, we ran the mechanical LHHM using our modi ed electrical model. We rst consider the pressure-volume (PV) loop. For the left ventricle, the loop plots pressure against volume throughout a cardiac cycle. The PV loop is used in studies to measure cardiac work output throughout the cycle (24).
Our modeled PV loop for the left ventricle is plotted against the PV loop from the control LHHM in Fig. 7. The comparison shows very little variation, ensuring negligible change to the mechanical contractionrelaxation cycle as a result of our modi ed Purkinje network. In particular, we nd the end diastolic volume (EDV), end systolic volume (ESV) and percentage ejection fraction (EF) compare well to literature as shown in Table 3 with the EF only slightly lower than the clinical average (43).
Moreover, peak left ventricle twist angle, the long axis (apex-to-base) shortening, and the mid-ventricular circumferential strain, all of which are markers of the contractile ability of the cardiac bers and in turn the ventricle's capacity to contract and effectively pump blood, were calculated for the control and for modi ed LHHM models (44). The results are shown in Table 4 which indicate full consistency between the two models. The long axis shortening and circumferential strain lie within the range of values from literature (45,46) while the left ventricle twist angle is slightly higher than the clinical average (44). According to the LHHM documentation the left ventricle twist angle is calculated as the apical rotation (47). Since most reports in literature have the twist angle as the apex-to-base rotation difference (44,48), the model results are expected have a slight overestimate.

Conclusion
Purkinje bers play a crucial role within the cardiac system to shape its ventricular activation sequence and corresponding ECG, as well as in synchronizing ventricular contraction and relaxation cycles. As such, they can no longer be pushed aside in in-silico models, as such models grow more and more clinically relevant. With the open availability of advanced algorithms that generate physiologically feasible networks, different Purkinje ber structures can be generated and adjusted to t and serve a wide range of cardiac models.
A state-of-the-art, electromechanical model like the LHHM is constantly changing, growing and adapting with the goal of providing insights needed to improve cardiac treatment and care. Having a healthy standard from which diseased models can be derived and compared is an obvious step in that direction. With the ECGs being a clinical standard in assessing cardiac function and investigating cardiac disorders, it was thus the aim of this study to produce pseudo-ECGs on the LHHM (state-of-the-art software) that match clinical normals, in preparation for ongoing developments and comparisons of diseased models and their ECGs. This we did by generating a two-stage Purkinje network with the rst stage stimulating endocardium sites of early activation prior to activating the remaining ventricular tissue. We then integrated our newly generated ber network into the electrical LHHM. The activation sequence produced from our modi ed electrical model was then successfully translated into pseudo-ECGs that were shown to be fully compatible with the average of clinical ECGs obtained from 39 healthy subjects whose information was made available on the PTB diagnostic database. This modi cation is thus an important extension to the capabilities of the out-of-the-box LHHM. Our subsequent analysis of the mechanical cardiac cycle with our modi ed model further con rmed its clinical consistency. The mechanical run of our modi ed model gives the same results as the LHHM (control) which means out model did not negatively affect the mechanical output but improved the electrical resemblance to clinical ECGs.
As a nal remark, adapting and evaluating the performance of our model under electromechanical disorders, such as left bundle branch block, right bundle branch block and dilated cardiomyopathy is of interest for future research. Moreover, to fully automate and optimize our Purkinje modi cation processes, for easy adoption by LHHM users, will also be the subject of our future research.

Declarations
Availability of data and material. All codes and materials used in this work have been referenced. All are freely available except the LHHM.
Competing interests. There are no competing interests to report in this manuscript.