In this study, we recruited six MI patients between 2018 and 2020 from Beijing Anshan Hospital and Dushuhu Hospital of Suzhou University. The study was approved by the Institutional Review Board of both hospitals. Cardiac magnetic resonance–late gadolinium enhancement (CMR–LGE) images of the patients were used to construct the heart models. A cardiac MRI was acquired using a 3.0 T scanner. Detailed image acquisition protocol has been published previously [30, 31, 32]. Table 1 provides detailed information on the images for each patient.

All analyses and measurements were performed using a custom software developed in MATLAB (Version: 2021a, Mathworks Inc., Natick, MA, USA). Two experienced experts manually segmented the epicardial and endocardial boundaries of the LGE images. The pixels between the boundaries were considered the myocardium. Then the modified Gaussian mixture model method was utilized to automatically identify the infarcted regions [33]. Finally, the full width at half max method was employed to further segment the infarcted tissue into the gray zone and core scar. The detailed process has been published previously [13, 33]. Following image segmentation, CardioViz3D [34] (INRIA, Sophia Antipolis, France) was used to interpolate the segmented low-resolution images to high-resolution images (approximately 0.4 mm). Then, the infarct tissue, including the core scar and gray zone, was interpolated using the log-odds method [35].

To quantitatively assess the impact of mesh resolution on the accuracy of electrical propagation in monodomain simulations, we introduced a cuboid of the same size used in the N-version benchmark [19] test. Table 2 provides detailed information about this cuboid, whose dimensions were 20 $\times$ 7 $\times$ 3 mm. Two different discretization methods were employed, which included uniform and adaptive tetrahedral meshes generated using Mesher in OpenCARP (opencarp.org) and 3-matic software (Materialize NV, Leuven, Belgium), respectively. As the adjustable parameters for different mesh resolutions varied between Mesher (shortest edge length) and 3-matic (maximal edge length), the mesh sizes for both discretization methods were different. The average mesh size of both meshes ranged from approximately 100 µm to 700 µm. Tables 3,4 (Ref. [19]) show the mesh statistics of both discretization methods.

First, we used the N-version benchmark to verify the accuracy of parameters used for both discretization methods. Then, we used a 10 mm cable and a slab with dimensions of 20 $\times$ 20 $\times$ 3 mm to examine the relationship between CV and mesh resolution. Table 2 provides the details of these three models.

To analyze the effect of varying model resolutions on the simulation results of patient-specific models, we generated six different models for each patient using Mimics (Materialize NV, Leuven, Belgium). Each model’s maximum and average edge length ranged from 0.45–1 mm and 0.3–0.65 mm, respectively. The upper limit of the model resolution was selected based on the previous benchmark studies of numerical convergence in cardiac electrophysiological modeling [19] and to test the effect of model resolution outside the theoretical convergence range on patient-specific ventricular model simulation results. Whereas its lower limit was selected based on the computer resources for model generation and high-preformation cluster.

After establishing the models, the fiber direction was specified in the mesh using a rule-based method [36]; The electrophysiological properties were assigned according to our previously published paper [13, 14]. Briefly, non-infarcted tissue was assigned to the human ventricular cell model published by ten Tusscher et al. [37]. Action potential remodeling in the gray zone, including changes in properties of fast delayed rectifier potassium, delayed rectifier potassium, L-type calcium, and sodium channels, was represented by reducing the corresponding currents (I${}_{\text{Kr}}$, I${}_{\text{Ks}}$, I${}_{\text{CaL}}$, and I${}_{\text{Na}}$) to 30%, 20%, 31%, and 38% of the normal values of the ten Tusscher ventricular myocyte model, based on previous experimental records [8, 38]. The core scar was modeled as an insulator that does not conduct electrical excitation. As shown in Table 2, the conductivities used for the cuboid were the same as in [19], but those for the cable, slab, and patient-specific model were set to match the clinically measured CV [10, 16, 39, 40, 41, 42].

The propagation of electrical activity in the heart model was simulated by solving a reaction-diffusion partial differential equation using FEM [43]. Electrical stimulations were performed using the openCARP simulation environment [44] on high-performance computers at the Dalian University of Technology, China, which was used to induce VTs as published previously [13]. All models were paced from 19 ventricular sites, including 17 sites on the left ventricle (LV), 1 near the right ventricular outflow tract, and 1 at the right ventricle apex, following the American Heart Association (AHA) Classification Standards [45]. After inducing reentry, arrhythmia was detected using a 10-second VT simulation based on the scheme proposed by Tong et al. [13].

The CV in the cable and slab simulations was calculated by taking the difference in activation times at locations 2.5 mm and 7.5 mm and 3.5 mm and 7.5 mm, respectively, and dividing it by the distance between the two points. CV values were calibrated by simulating wavefront propagation following stimulation at the center of a slab model (20 mm $\times$ 20 mm $\times$ 3 mm) with longitudinal fiber orientations along the X-axis. To achieve the same CV in models with different mesh resolutions, conductivities in all models were adjusted to match the CV of the highest-resolution model.

The VT locations in the clinic were measured as part of the electrophysiological study performed during implantable cardioverter defibrillator (ICD) implantation or estimated by a senior electrophysiologist via electrocardiogram (ECG). To compare the locations quantitatively via clinical assessment and simulation, the VT locations were assigned to the 17 segments following the American Heart Association Classification [45]. If the segment of the predicted VT location in the simulation was the same as the clinical measurement, it was considered clinically relevant reentry, or otherwise it was treated as a clinically irrelevant VT.

Fig. 1C displays the activation times at point P8 for solutions with different time steps for adaptive and uniform meshes of the cuboid (Fig. 1A,B). The activation time at point P8 with a uniform mesh had below 10% error at a spatial scale of approximately 340 µm (green line in Fig. 1). When a spatial scale of uniform mesh size of 600 µm (718.2 $\pm$ 124.1 µm) was used (Table 3), the error was 65.1% compared to the assumed converged solution of 43 ms as reported in reference [19]. The activation time displayed nonmonotonic behavior when the mesh resolution was increased to 90 µm (108.5 $\pm$ 18.6 µm).

Fig. 1.

Model discretization with adaptive (A) and uniform tetrahedral (B) mesh. Activation times at point P8 for solutions with dt = 25 µs and 10 µs for adaptive meshes, and dt = 25 µs for uniform meshes (C). The ground truth of 43 ms is the assumed converged solution, as reported in reference [19].

For simulations with adaptive meshes (red and blue lines in Fig. 1C), we observed minor differences between time steps of 25 and 10 µs. The error in the relative difference in activation time at P8 with adaptive tetrahedral meshes was within 10% compared to the assumed converged solution of 43 ms [19], although the simulation results with adaptive meshes also displayed nonmonotonic behavior when the mesh resolution was changed from an average value of 700 µm to 100 µm.

Table 5 summarizes the volumes of normal myocardium, gray zone, and core scar of the reconstructed cardiac models. The average ventricular volume was 174.3 cm${}^3$. The percentage of myocardial volume in non-infarcted areas and the gray zone was approximately 88.5%–96.9% and 2.5%–8.7%, respectively, with the core scar comprising 0.5%–4.9%. The variation in volume for each tissue (non-infarct tissue, gray zone (GZ), and core scar) in the models with different resolutions was within 0.01 cm${}^3$.

Table 6 presents the statistical characteristics of the generated tetrahedral models with different mesh resolutions for all six patients. Res1–Res6 in the first column represent different model resolutions; the second column displays the value of the maximal edge length of the tetrahedral mesh for each model in the Mimics software; the third and fourth columns show the average number of vertices and elements for all six tetrahedral models with the corresponding resolution, respectively; the fifth column represents the average edge lengths for all six models with different maximum edge lengths; and the last column shows the simulation time for each time step of 10 ms.

For the cable simulation, the results display nonmonotonic behavior similar to that of the cuboid (Fig. 2A). The CV in healthy tissue was the fastest (47.19 cm/s) with a uniform mesh size of 150 µm (173.8 $\pm$ 30.2 µm). The maximal error of CV was 23.4% compared with the value in the finest uniform mesh size of 100 µm (115.9 $\pm$ 20.2), and the CV decreased linearly between the uniform average mesh size of 300 µm and 700 µm (red line). Similar behavior was observed in the GZ (blue line). The results of the slab simulation exhibited similar nonlinear behavior compared to the cable simulation. However, the absolute longitudinal CV decreased for healthy tissue and the GZ (Fig. 2B). The absolute and relative differences in CV were 5.1 cm/s and 12.4%, respectively, between average mesh sizes of 142 and 698 µm (red circle in Fig. 2) in the slab simulation.

Fig. 2.

CV variation at different mesh resolutions. (A) CV of healthy tissue and gray zone in the cable model. (B) Longitudinal and transverse CV of healthy tissue and gray zone in the slab model. CV, conduction velocity.

Table 7 summarizes the longitudinal CV measured in the slab meshes with different resolutions. Res6 has the highest CV of 41 cm/s, whilst Res1 has the slowest CV of 35.9 cm/s for the non-infarct tissue, making the difference 5.1 cm/s. For the gray zone, the difference in CV between the two-resolution models was 5.7 cm/s, which was higher than that between the non-infarct tissues. To compensate for the lower CV caused by the coarsening of mesh resolution, the CV in models with Res1–Res5 was adjusted to match the highest-resolution CV (Table 7).

Table 8 summarizes the overall simulation results and the location of the reentries induced in models with different mesh resolutions for all six patients. The VT inducibility and induced VT morphology and location were very similar in the models with resolutions from Res3 to Res6. Using Res6 as a criterion, the results in Res1 and Res2 models varied considerably.

The simulation results in models with different resolutions were much more consistent for patients without clinically observed VT. No VT was induced at any sites in the models PAT02, PAT04, and PAT05 for all mesh resolutions before and after CV adjustment. For PAT01, one pacing site in the models with Res4 and Res6 induced reentry. However, the VT-inducing ratio (5.3%) was very low, suggesting that the patient had a low probability of inducing VT. For patients with clinically observed VT, the simulation results in models with different mesh resolutions were more diverse than those without clinically observed VT (Table 8).

In PAT03, we monitored the VT events four times in the clinic and found that the VT was located at the middle lateral wall of the LV (Segment 7 in AHA 17 segments). In the Res1 and Res2 models, three pacing sites induced a type of VT located at the posterior wall (Fig. 3A), which was unrelated to the location measured in the clinic. Though the VT induction ratio is comparable to other higher resolution models, no pacing site-induced VT related to clinical location (Fig. 3B). This VT morphology was also present in the Res3 and Res4 models (Fig. 3C,D, respectively), although at significantly smaller proportions (33.3% and 28.6%, respectively).

Fig. 3.

VT morphology located at Segment 7 (clinically irrelevant) induced in simulation for PAT03 with different mesh resolutions. VT induced in (A) Res1, (B) Res2, (C) Res3, and (D) Res4 models. The red arrowhead indicates the direction of VT propagation. The color scales in (A–D) indicate activation times, and the black areas represent the core scar (without any electrical propagation). VT, ventricular tachycardia; CL, cycle length.

The results of the Res3–Res5 models are more consistent with the Res6 model (Table 8). They all induced clinically relevant VT (Fig. 4), with Res4 and Res6 inducing a maximum of five sites. Res6 did not induce any VT unrelated to the location measured in the clinic. Although Res3–Res5 had VTs unrelated to the clinical measurement, the percentage of incorrectly predicted VTs was low (20%–33.3%). When the CV was adjusted for models, the incorrectly predicted VTs were reduced in almost all models (Table 9). After the CV adjustment, Res3 was the only model with decreased accuracy. Though the incorrectly predicted VTs in segment 10 were reduced, two new incorrectly predicted VTs emerged. For the Res4 model, the number of clinical-related VTs increased from five to eight (Table 9). Although a new VT was induced at a certain pacing site that did not exist in the model with original conductivity, it cannot be considered clinically relevant owing to its low occurrence. For the Res5 model, the number of clinical-related VTs increased, while the VT unrelated to the clinical measurement disappeared. In general, the accuracy of VT prediction in the Res1, Res2, and Res3 models was not high, even after CV adjusting. Fig. 4 shows the VT morphology located at segment 7 for PAT03 with different mesh resolutions.

Fig. 4.

VT located at Segment 7 (clinically relevant) induced in simulation for PAT03 with different mesh resolutions. VT was induced in the (A) Res3, (B) Res4, (C) Res5, and (D) Res6 models. The red arrow indicates the direction of VT propagation, the color scales in (A–D) indicate activation times, and the black areas represent the core scar (without any electrical propagation). VT, ventricular tachycardia. CL, cycle length.

For PAT06, although all models with different resolutions induced clinically relevant VT (Fig. 5), the percentage of pacing sites that included VT was low for the Res1 and Res2 models (5.3% and 15.8%, respectively). However, for the Res3–Res5 models, the number of pacing sites that induced VT was very close to the highest resolution ($\pm$ pacing sites). After adjusting the CV, all models simulated the same VT morphology at the same clinically relevant location. The Res2–Res5 models predicted a decrease in the number of VTs, while Res1 predicted an increase in the same (Table 9).

Fig. 5.

VT morphology located at segment 7 (Clinically relevant) induced in simulation for PAT06 with different mesh resolution. VT induced in (A) Res1 model, (B) Res2 model, (C) Res3 model, (D) Res4 model, (E) Res5 model, and (F) Res6 model. The red arrow indicates the direction of VT propagation, the color scales in (A–F) indicate activation times and the black areas represent core scar—there is no electrical propagation there. VT, ventricular tachycardia; CL, cycle length.

A limitation of this study is the small sample size, consisting of only six patients. This is due to the challenges in obtaining good-quality MRI images and patient follow-up data. Despite the small patient dataset, we believe the conclusion remains unaffected. Another limitation is that the smallest mesh size is approximately 300 µm, a constraint imposed by the 64GB memory capacity of the computer utilized for model generation. As the simulation results in cable and patient-specific models with Res5 and Res6 exhibited no significant differences, we believe that models with resolutions higher than 300 µm will not be adversely affected.