Edited by: Javier Saiz, Universitat Politècnica de València, Spain
Reviewed by: Olivier Bernus, Université de Bordeaux, France; Bradley John Roth, Oakland University, United States
This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Non-ischemic dilated cardiomyopathy (NIDCM) is characterized by enlarged ventricular cavity size, and impaired ventricular systolic function, which is not a consequence of myocardial ischaemia. Patients who present with NIDCM are known to be at high risk of sudden cardiac death (SCD), with an estimated 20% mortality rate over 5 years (Gulati et al.,
Further risk stratification in NIDCM based on LGE-CMR is therefore vital, yet is challenged by the lack of knowledge of both the underlying fibrosis micro-structure, as well as its implications for arrhythmogenic mechanisms. This is partially due to the limitations of current clinical scan resolutions, which are an order of magnitude larger than that required to resolve micro-structural fibrosis (Schelbert et al.,
In-light of the limitations of LGE in directly quantifying fibrosis density in non-ischemic disease, it may be possible to enhance risk stratification in NIDCM patients by combining LGE assessment along with functional electrical measurements, such as 12-lead ECG data or invasive electroanatomical mapping recordings (Betensky et al.,
Computational modeling provides a useful tool to better understand the mechanistic role of cardiac micro-structural remodeling in the initiation and maintenance of ventricular arrhythmias in the context of mid-wall fibrosis. Recently, detailed computational models have been applied to understand the arrhythmogenic interaction of electrical wavefronts and micro-fibrotic structures in a variety of pathological conditions such as modeling of atrial fibrillation (Roney et al.,
An anonymized LGE-CMR image set of a patient presenting with NIDCM and LGE was acquired from the Royal Brompton Hospital, using a previously described protocol (Gulati et al.,
Computational model creation from an LGE-CMR image.
In Figure
Interstitial and replacement fibrosis were modeled by modifying mesh edges and triangles, respectively, with a single parameter α controlling the density of fibrosis in LGE. The parameter α was selected from the range 0−1 in steps of 0.1, which yielded models with varying levels of fibrosis. Additionally, each model's level of fibrosis was related to the normalized image intensity
with
For replacement fibrosis, the probability of a mesh triangle being fibrotic was
where α is the global fibrosis density. Fibrotic triangles were removed from the mesh, making them electrically inert and non-conducting.
Interstitial fibrosis was represented by a network of random fibrotic clefts (Costa et al.,
where θ is the angle between the triangle edge and the local myocardial fiber direction. The cos4(θ) term in the equation greatly decreases the likelihood of fibrosis not aligned with the myocardial fiber direction. This ensures a highly anisotropic fibrosis distribution, which is what we would expect in the case of interstitial fibrosis.
Model tissues within circuits of interstitial fibrosis were removed, since they were electrically isolated. Example models with both interstitial and replacement fibrosis at various densities are shown in Figure
Simulated fibrosis micro-structures with varying type and density. Black lines are separating edges through which current cannot cross, whereas white areas are electrically inactive and non-conducting.
Conductivity values in non-enhanced areas were tuned to match conduction velocities (CV) from NIDCM (Anderson et al.,
Both normal conduction (NC) and reduced conduction (RC) were considered in LGE. In the case of NC the same conductivity values were used as in non-enhanced areas. In the case of RC, model LGE areas were assigned one of four reduced CVs based on image intensity. Regions in the intensity range 0–25 and 25–50% above threshold had transverse conductivity (CVT) reduced by 25 and 50%, respectively, with normal fiber direction conductivity (CVF). Finally, regions in the intensity ranges 50–75 and 75–100% above threshold had CVF reduced by 25 and 50%, respectively, and CVT reduced by 50%. These four reduced sets of CV values are consistent with (Anderson et al.,
We constructed a set of model fibrosis micro-structures, consisting of variations in type (replacement, interstitial), conductivity (normal, reduced) and 10 levels of maximum density α (0.1–1.0). Ten random realizations of fibrosis were made for each combination of conductivity, fibrosis type, and fibrosis density, along with a control model with normal conductivity and no fibrosis. This gave a total of 401 models (Figure
Electrical activity was simulated by the standard monodomain formulation, with ionic currents represented by the Ten Tusscher and Panfilov (
Models were paced primarily from a site on the LV septal endocardium, approximately halfway down the extent of the LGE (see Figure
Pseudo-electrograms were calculated by estimating the extracellular potential (Bishop and Plank,
where Ω is the cardiac domain, σ
Activation patterns and transmural activation times were generated with a fixed stimulation protocol, consisting of 5 beats with coupling intervals (CIs) 500, 340, 250, 210 ms.
Reentry inducibility was tested with a dynamic algorithm, with up to six stimuli with variable CI. The algorithm consists of the components: reentry detection, coupling interval selection, and stimulus capture detection.
For each stimulus 600 ms of electrical activity were simulated, and electrical activations were recorded. Reentry was detected if any activations were present after 170 ms in a region of interest (ROI). For all experiments the ROI was defined by a circle with radius 3 mm, centered at the stimulus site, and excluded fibrotic areas.
The coupling intervals (CI) between stimuli were determined algorithmically for every stimulus after the first. First the CI of 200 ms was attempted, and if the resulting stimulus generated a new wave, then this CI was kept. Otherwise the effective refractory period (ERP) at the stimulus site was determined, and the CI was set to ERP + 1 ms. This ensured that electrical stimuli were delivered rapidly enough to destabilize the model, and yet slowly enough that each stimulus generated a distinct wave. An effective lower bound of 200 ms was placed on the CI, as pacing rates faster than this can generate reentries under non-pathological conditions (Cao et al.,
Knowledge of whether a stimulus generated a new activation wave was required for the coupling interval selection. This was accomplished by tracking a wave of activations for 20 ms after a stimulus was applied, with a stimulus determined as capturing if activations could be tracked up to 20 ms. Initially, all activations inside the reentry detection ROI were found 1 ms after stimulus delivery. For each subsequent millisecond new activations were located inside a set of tracking ROI, which consisted of model tissue that was within the distance that an activation wave could travel in 1 ms, assuming the maximum CV of 84 cm/s.
Reentries were classified by visual inspection of videos of the simulated transmembrane potential. A reentry was classified as a rotor if an organizing center was present, consisting of nearly simultaneous activation and repolarization, as micro-reentry if a narrow conducting channel was involved, and as macro-reentry if activation followed a large circuit. When multiple reentry types were present in a simulation, the reentry type of the earliest circuit was used. Examples of the three reentry types are shown in Figure
Examples of three canonical reentry mechanisms.
Figure
Incidence of reentry mechanisms by fibrosis texture, density, and with normal and reduced conductivity in the LGE.
The incidence of each reentry mechanism is shown in Figure
We examined the effects of fibrosis type and pacing rate on local activation patterns around the region of LGE. Zoomed in views of the activation maps at CIs, 500, 340, 250, 210 ms, for fibrosis densities 0.5, 0.8, and 1.0, along with the control case, are displayed in Figure
Activation maps around the LGE for three types of fibrosis micro-structure. Red lines are isochrones spaced at 10 ms intervals. The stimulus site is shown in green in the fibrosis-free control model (bottom row).
We tested the ability of TAT to predict reentry inducibility as a surrogate for knowledge of the fibrosis micro-structure. In particular we considered three metrics, TAT500, TAT210, and ΔTAT. TAT500 and TAT210 refer to TAT at CIs 500 and 210 ms, respectively, whereas ΔTAT refers to the difference between TAT210 and TAT500. Models for which reentry was inducible had significant higher mean values of the three metrics (Mann–Whitney
Using cut-off values of the TAT metrics, we predicted reentry inducibility using all of the model micro-structures. The results are summarized in Figure
ROC curves for predicting reentry inducibility based on transmural activation times (TAT). TAT500 and TAT210 are the transmural activation times with 500 and 210 ms coupling intervals respectively, while ΔTAT is the difference. The green dot shows the point of maximum sensitivity and specificity. AUC stands for area under curve.
We examined the effects of pacing and measuring TAT from alternate locations (see Figure
Effect of pacing location on transmural activation times and reentry incidence. Transmural activation times are given as a mean (solid line) and standard deviation (transparent area) for four coupling intervals (CI), taken over 10 random realizations of each fibrosis density. The reentry incidence is overlaid in red for comparison. All results were obtained with replacement fibrosis texture and reduced conductivity. The results for site 1 are the same as in Figure
Our study has highlighted the role of fibrosis density and local conduction slowing for the induction of reentry in NIDCM. Furthermore, we identified distinct reentry mechanisms related to the underlying fibrotic substrate, and found that variations in micro-structure density and conductivity had a measurable effect on activation patterns. This allowed for the accurate prediction of reentry inducibility in a wide variety of fibrotic micro-structures with pacing from the endocardium halfway through the extent of the scar.
To the best of our knowledge, our paper is the first to develop a computational model of fibrotic scarring in NIDCM. The main feature of our model is the incorporation of random fibrosis micro-structures. These micro-structures cannot be imaged with current clinical modalities, and were therefore approximated by our computational model.
Micro-structures have explained reentry formation in theoretical studies (Alonso and Bär,
We modeled reduced passive tissue conductivities within fibrotic areas, based on evidence of reduce gap junction protein expression in the context of fibrosis and dilated cardiomyopathy (Kostin et al.,
We modeled the excitability of fibrotic tissue in NIDCM as normal, which is supported by experimental studies (Anderson et al.,
We did not include any models of myocyte-fibroblast coupling (MacCannell et al.,
Finally, we did not include the effects of APD prolongation and altered transmural heterogeneity as were observed in Glukhov et al. (
Previous studies have shown that fibrosis in NIDCM plays a role in both the initiation of block, and in the slowing of propagation (Anderson et al.,
By combining the known consequences of fibrosis in NIDCM into a computational model, we were able to simulate three reentry mechanisms, and analyse the relative contributions of fibrosis micro-structure and RC. We observed that macro-reentries involved large areas of functional and anatomical conduction block, which were more likely to form when the fibrosis was sufficiently dense, especially for replacement fibrosis with NC and interstitial fibrosis with RC.
We observed that rotor reentries formed when relatively large wavefronts pivoted and folded back onto an area of transient conduction block. The anisotropic properties of interstitial fibrosis facilitated this. Relatively free conduction in the transverse direction allowed for the preservation of large wavefronts in the fibrotic areas, whereas restricted conduction in the fiber direction caused zones of transient block which could be reentered by the pivoting waves. Indeed, most of the rotor reentries occurred with interstitial fibrosis. RC also played a key role in rotor formation, as only as single rotor was observed with NC.
We observed micro-reentries when propagation was constrained to a narrow pathway, and yet was not completely blocked. This occurred predominantly with replacement fibrosis and RC for the density range 0.6–0.9, at which micro-reentries were more common than macro-reentries. We did not observe micro-reentries with interstitial fibrosis. This was most likely due to the relatively unhindered lateral conduction of interstitial fibrosis which did not support the formation of narrow channels required for micro-reentry.
Previous simulation studies (Alonso and Bär,
All reentries were initiated in the presence and absence of RC. The possibility that reentry in NIDCM is possible due to the structural effects of fibrosis alone, without RC, is also supported by previous simulation studies (Kazbanov et al.,
A comparison of our simulated TAT values to literature provides a preliminary validation of our NIDCM fibrosis model. In particular, for our control case with no fibrosis and NC, we observed a TAT of 25 ms, which is within the range 26 ± 14 ms observed in healthy controls by Vassallo et al. (
We found that activation in the fibrotic region followed one of two patterns, direct or compartmentalized. These same patterns were observed by Betensky et al. (
Our ROC analysis indicates that TATs are predictive of reentry. This is supported by Betensky et al. (
We considered alternative pacing and TAT recording locations and measured the reentry incidence and TAT values at these sites. At the three sites in the center of the scar, the reentry incidence was similar and the TAT values correlated with scar density. At the off-center locations the reentry incidence was lower, with the site furthest away from the center of the scar experiencing the lowest incidence. This is expected since waves approaching the scar from the off-center locations traveled in a vector that was more greatly aligned with the myocardial fiber direction than for waves traveling from the central pacing locations. This meant that the effective CV of waves originating form the off-center locations was higher, and conduction block less likely.
Furthermore, TAT was only weakly correlated with the scar density at the off-center pacing sites. This was due to the effects of waves traveling around the scar, which were more likely to be the first to reach the TAT recording locations in the off-center setups. Since these waves traveled around the scar rather than through it, they experienced less fibrosis induced delay. For the pacing cite furthest away from the center of the scar, the scar transmurality was also substantially less, which lead to less delay for waves traveling through the scar than from the other pacing sites.
The DANISH trial (Køber et al.,
We demonstrated good performance of ΔTAT for predicting simulated reentry induciblity for a single scar macrostructure with a variety of miscrostructures. Further experiments with a variety of macrostructures are warranted. Our results with alternative pacing and TAT recording locations suggest that finding a location central to the scar is crucial for inferring the potential for reentry from TAT measurements. Furthermore, the method is robust to small changes in the pacing/recording locations. This is important for any potential clinical application, where less precision in determining locations can be expected.
Another potential data source complementary to LGE-CMR is T1 mapping, which allows for absolute comparisons to be made, and has been shown to be an independent predictor of arrhythmia (Chen et al.,
Computational models of NIDCM, building on the methodology presented here, have the potential to predict SCD in a manner similar to that which has recently been demonstrated for ischemic disease (Arevalo et al.,
All models were based on a single image. Consequentially our results do not account for geometric variability, such as wall thickness, scar transmurality, and scar extent. These factors may effect both TAT and reentry inducibility, and therefore the ability to infer reentry susceptibility from TAT measurements. Furthermore, our simulated TAT differences were a consequence of electrical restitution properties, which were governed by the Ten Tusscher and Panfilov (
As shown in section 3.5, the choice of pacing location is critical when making TAT measurements. In the clinical study of Betensky et al. (
Simulations were performed in 2-D, which neglected the 3-D structure of ventricular tissue, myofibers, and interstitial clefts. As a consequence we restricted our study to the generation of reentry and not the sustenance, in which additional electrical pathways in 3-D may be important. Furthermore, micro-structural effects on TAT values could be less pronounced in 3-D at comparable fibrosis levels. This is because the presence of a 3rd dimension would likely increase the degree of connectivity among different paths through the scar, leading to quicker fill-in by multiple propagating activation wavefronts and hence more stable conduction. At the same time source-sink mismatches could be more pronounced in 3D due to the extra electrotonic load on a point in 3-D tissue. This could make conduction block and hence reentries more likely to occur.
Furthermore, the modeled 2-D fiber architecture did not account for the effects of transmural variation in fiber helix angle. Including this effect into a projection of the fibers onto a 2-D short-axis plane would cause the effective conductivity in the in-plane fiber direction to be less at the epicardial and endocardial edges. We expect that including this in our models would act to shorten the electrical wavelength around the scar, potentially increasing the reentry incidence.
Only two levels of conductivity reduction were tested: none and the heterogeneous scheme described in section 2.4. This was due to the lack of experimental data to inform the conductivity values, and to limit the number of model combinations.
Finally heterogeneous blocks of fibrosis corresponding to image pixels were present in our models, which have the potential to affect the probability of reentry (Kazbanov et al.,
All authors contributed to the interpretation of the data, drafting, and revising the manuscript, and approved the final version of the manuscript. The original study design was made by GB, MB, SP, BH, and discussed with the other authors. GB ran the simulations and analyzed the resulting data.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors thank Hermenegild Arevalo for his advice with regards to the design of pacing protocols.
The Supplementary Material for this article can be found online at: