Study populationThis study was approved by the Institutional Review Board of the Severance Cardiovascular Hospital, Yonsei University Health System, and conducted in accordance with the principles of the Declaration of Helsinki. All patients in the Yonsei AF Ablation Cohort Database (ClinicalTrials.gov Identifier: NCT02138695) provided written informed consent for their clinical data to be used in digital twin studies. We followed the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guideline42.To select patients with AF for whom CT and EAM could be performed, this study was conducted on those who underwent AMD treatment for recurrent or symptomatic non-sustained AT within 3 months after AFCA. Between April 2014 and October 2021, 681 patients who underwent AFCA and were eligible for a digital twin study were initially enrolled. The final cohort comprised individuals who underwent AMD treatment within 3 months and were followed up for up to 1 year from the commencement of AMD therapy. Individuals who abstained from AAD treatment post-AFCA, patients for whom AAD was initiated more than 3 months after ablation, patients treated with drugs other than AMD, patients with a follow-up duration <1 year, and cases without additional analysis in the digital twin study were excluded. Finally, 115 patients were included in the analyses (Fig. 2).Clinical substrate mapping and AF ablationAt AFCA initiation, we created LA substrate maps (comprising bipolar voltage and local activation maps) using the EnSite NavX system, obtained with a multielectrode catheter (AFocus, Abbott, Chicago, IL, USA) during high right atrial (RA) pacing, with a cycle length of 500 ms. If AF persisted, we performed internal cardioversion using biphasic shock (2–20 J) with R-wave synchronization (Lifepak12, Physiocontrol Ltd., Redmond, WA, USA) to restore sinus rhythm. Cases in which sinus rhythm could not be restored, PVI was performed, followed by mapping after sinus rhythm restoration. We collected bipolar electrogram data from the LA surface and obtained more than 500 points per patient.An open-irrigated tip catheter (FlexAbility, Abbott Inc.; Coolflex, Abbott Inc., Minnetonka, MN, USA; 30–35 W; 47 °C; TactiCath, Abbott Inc.; and ThermoCool SmartTouch, Biosense Webster Inc.) was employed for AFCA. The ablation endpoint at each site was defined as an average impedance drop >10% of the baseline or a >80% decrease in the local electrogram voltage amplitude. We achieved circumferential PVI (CPVI) using a bidirectional block in all patients. Additional lesion formation beyond CPVI was determined at the operator’s discretion.Follow-up strategy and AMD prescriptionAfter the procedure, all patients were monitored regularly following the same schedule as that followed for previous patients who underwent AF ablation in our center43. AMD was prescribed, if no contraindications were present, to patients with early recurrence of AF or AT within 3 months, symptoms of suspected arrhythmia, or ineffective cardioversion during the procedure. Typically, AMD was administered at a dose of 200 mg once daily for 2–6 weeks, followed by a reduction to 100 mg if rhythm control was achieved, which was then maintained. Patients who underwent post-AFCA AMD treatment were followed up at 1, 3, and 6 months after the procedure and every 6 months thereafter.Creating a digital twin of the LAWe created a digital twin of the LA by merging pre-ablation CT images with EAM data obtained during AFCA.3D Geometry formationThe structural components of the geometry were represented as triangular meshes, with each node forming a human atrial myocyte (Fig. 1b). First, triangular meshes were developed based on individual LA CT images. These 3D mesh surfaces consisted of ~400,000–500,000 nodes, with a mean spacing of 235.1 ± 32.1 μm between adjacent nodes. Throughout the AFCA procedure, we collected bipolar electrogram data from >500 points on the LA surface using a circular mapping catheter and an Ensite NavX system (Abbott Inc., Chicago, IL, USA), during a paced rhythm with a cycle length of 500 ms. We aligned the coordinates of the EAM with individual CT images, producing a clinical EAM. Additional details regarding the procedure with aligning EAM with CT images were described in the previous study44.Human atrial myocyte model and action potential propagationWe used a modified Courtemanche-Ramirez-Nattel (mCRN) model45,46 to characterize the system, which mathematically represents the various ion currents within human atrial myocytes, as shown in Eq. (1).$$\begin{array}{l}{I}_{{ion}}={I}_{{Na}}+{I}_{{CaL}}+{I}_{{to}}+{I}_{{Kur}}+{I}_{{Kr}}+{I}_{{Ks}}+{I}_{K1}+{I}_{{K}_{{Ach}}}\\\qquad\;\; +\,{I}_{{NaCa}}+{I}_{{NaK}}+{I}_{b,{Na}}+{I}_{b,{Ca}}+{I}_{p,{Ca}}\end{array}$$
(1)
Where Iion denotes the ionic current of the atrial myocyte, comprised of various individual ion channel currents—including sodium (INa), potassium (Ito, IKur, IKr, IKs, IK1, IKAch), calcium (ICaL), exchanger (INaCa), pump (INaK, Ip,Ca), and background currents (Ib,Na, Ib,Ca). This model served as the foundation for generating action potentials and simulating wave propagation in atrial myocytes.Assuming each cardiac cell represents a single node, a triangular array was created with nodes representing cell-to-cell connections. To computationally model cardiac action potential propagation through the atrial wall, we used the following reaction-diffusion equation47 referred to as Eq. (2).$$\frac{\partial {V}_{m}}{\partial {\rm{t}}}=\nabla \,\cdot\, D\nabla {V}_{m}-\frac{{I}_{{ion}}+{I}_{{stim}}}{{C}_{m}}$$
(2)
Where Vm (V) denotes the membrane potential, Cm (F/m²) represents the membrane capacitance per unit area, D (m2/s) is the diffusion coefficient, and Istim (ampere/meter²) refers to the stimulation current. To simulate the reaction-diffusion system in 3D, we used models constructed along a generalized finite difference scheme.Electrophysiological and histological characterizationAfter geometry formation, the electrophysiological characterization of each node involved determining parameters, such as voltage, fibrosis state, fiber orientation, and conductivity (Fig. 1b). Voltage interpolation was performed for clinical voltage data using the inverse distance weighting method as shown in Eq. (3)48$${W}_{{ij}}=\frac{{{d}_{{ij}}}^{-a}}{{\sum }_{k}^{{n}_{j}}{d}_{{kj}}},{R}_{j}=\mathop{\sum }\limits_{i=1}^{{n}_{j}}{w}_{{ij}}{R}_{{ij}}$$
(3)
Where W represents the weighted average of neighboring values, i and j indicate the known and unknown values of the points, dij− a is the distance between known and unknown points, Rij represents the value of the known point, and Rj indicates the interpolated value at the unknown point j. The interpolation process produced the virtual voltage data with an amplitude within a 10-mm radius from the region of interest.The fibrotic state was determined by applying the obtained virtual voltage values to a probability function, to distinguish fibrotic cells from normal cells. To determine fibrosis status (yes/no) for each node, we used the following nonlinear equation (Eq. (4)) between the bipolar voltage and probability of fibrosis:$${P}_{{fibrosis}}=\left\{\begin{array}{cc}1 & V\, <\, 0\\ \frac{1}{100}* (-40{V}^{3}+155{V}^{2}-206V+99.8) & 0\le V\le 1.74\\ 0 & 1.74\, <\, X\end{array}\right.$$
(4)
where Pfibrosis is the probability that there is fibrosis at a given node, and virtual V is the bipolar voltage at that node within the 0–1.74 mV range. If V is >1.74 mV, Pfibrosis converges to zero. This was developed by comparing the predicted percentage of fibrosis across the 3D atrial model with pre- and post-ablation fibrosis data. For each node, the probability of fibrosis calculated based on the clinically acquired bipolar voltage data was compared against a random number between 0–1. In fibrotic cells, ion currents—including the inward rectifier potassium current (IK1), L-type calcium current (ICaL), and sodium current (INa)—were reduced by 50%, 50%, and 40%, respectively12.For fiber orientation, we used an atlas49 based on the fiber orientation of the LA surface and compared it with the CUVIA program to adjust for patient-specific geometries. Representative vectors were drawn, and fiber orientations at surrounding nodes were determined through interpolation around these representative vectors. A vector aligned with the myocardial fiber direction was generated at each point in the heart. Conductivity was set based on orientation, differentiation between the longitudinal and transverse directions, and fibrotic and non-fibrotic tissues12. We defined the longitudinal conduction velocity as that in the same direction as the vector and the transversal conduction velocity as that in the perpendicular direction to the vector. The conductivity of the model was applied at 0.1264 S/m (non-fibrotic longitudinal cell), 0.0546 S/m (fibrotic longitudinal cell), 0.0252 S/m (non-fibrotic transverse cell), and 0.0068 S/m (fibrotic transverse cell).Finally, to synchronize the clinical and virtual LAT maps and determine the conduction velocity, we displayed the clinical LAT map on the CUVIA program screen. An experienced investigator synchronized the virtual LAT map to match its appearance to that of the clinical map. To assess inter-observer variability, we compared LATs synchronized by two different investigators across 10 cases and fiber orientation maps created by three different investigators. After inducing virtual AF, we compared the outcomes. We found that all 10 cases resulted in the same final rhythm (AF, AT, or termination), regardless of the synchronization or fiber orientation method used. Fleiss’ Kappa coefficient for inter-observer variability was calculated to be 1, indicating perfect agreement. By adjusting the diffusion coefficients of the 3D model to match the conduction velocity of the virtual LAT map with that of the clinical LAT map, we were able to induce virtual AF and interventions in the digital twin.Virtual PVI and AMD interventionUsing the CUVIA digital twin, circular lesions of 2 mm width were created on both sides of the PVs (Fig. 1c). We conducted PVI at the antral levels. In the presence of AF, which is characterized by distinct ion current adaptations compared with sinus rhythm, adjustments were made to the conductance or concentration of specific ion channels (Supplementary Table 3). Considering the increased recurrence risk (early recurrence, symptoms, or challenges in achieving sinus rhythm during the procedure) among enrolled patients, we tailored the diffusion coefficient of the control model to ensure the persistence and maintenance of AF in all simulations during the observation period (32 seconds).To conceptualize the conditions under which AMD acts within the body at subtoxic ranges, we defined low, high, and toxic doses of AMD as 1.6 μM (minimal effective concentration), 3.9 μM (maximal effective concentration), and 8.0 μM (toxic concentration) respectively (Fig. 1d), based on the therapeutic range of AMD50. Given the role of AMD as a multiple ion channel blocker, we investigated the functional blockade of ion channel conductance at these concentrations. The degree of functional blockade was evaluated using Hill’s equation8,33, expressed as Eq. (5).$$\theta ={[1+{({{IC}}_{50}/D)}^{{nH}}]}^{-1}$$
(5)
Where θ denotes the degree of channel blockade (ranging from 0 to 1), IC50 represents the half-maximal inhibitory concentration, D is the free drug concentration, and nH is the Hill coefficient. The ion channel specific IC50, Hill’s coefficient, and corresponding references, along with the ion current settings for low (1.6 μM), high (3.9 μM), and toxic (8.0 μM) AMD concentrations, are summarized in Supplementary Table 3.Evaluation of electrophysiologic parameters depending on virtual AMDTo assess electrophysiological changes under AMD treatment, we conducted pacing at the earliest activation site (EAS) with a cycle length of 500 ms (Fig. 1e). We measured APD90 and peak upstroke velocity across baseline, low-, high-, and toxic-dose AMD scenarios near the activation site. Subsequently, we compared the mean values and assessed whether a dose-dependent trend was observed.Next, to assess the effectiveness of AMD in patients, all conditions (baseline, low-, high-, and toxic-dose AMD) were subjected to ramp pacing to induce AF (Fig. 1g). To induce virtual AF, we ramped pacing around the EAS for 11.52 s with 8 beats/cycle. The pacing started at 200 ms and was decreased at 10 ms intervals until it reached 120 ms.We created the Smax and DF maps for the aforementioned scenarios as shown in Fig. 1g. To determine the Smax, the APD90 and DI were measured at each node during a pacing period of up to 3 beats after a wave break during rapid pacing (200 ~ 120 ms)26. Smax was calculated as the maximum slope of the APD90 restitution curve and defined for all nodes (>400,000) in the LA model. The nonlinear fitting of the APD90 and DI was calculated using the following correlation equation, Eqs. (6) and (7).$$y\left({Action\; potential\; duation}\right)={y}_{0}+{A}_{1}\left(1-{e}^{-\frac{{DI}}{{\tau }_{1}}}\right)$$
(6)
$${slope}=\left(\frac{{A}_{1}}{{\tau }_{1}}\right)* {Exp}\left(-\frac{{DI}}{{\tau }_{1}}\right)$$
(7)
Where y0 and A1 are free-fitting variables, and τ1 is a time constant. In each patient, we obtained the Smax value of each node.The DF was derived from a Fourier transform of the action potentials for 6 s (17–23 s after the initiation of ramp pacing) at each node15. We evaluated the DFs for all nodes of the LA model.After virtual PVI, the remaining LA was divided into 6 regions (septum, anterior wall, left atrial appendage, left lateral isthmus, posterior wall, posterior inferior wall) per patient to perform regional analysis (Fig. 1f)25. In each region, we analyzed the average Smax and DF values of each regions (Fig. 1g). We also analyzed the highest and lowest mean Smax and DF values among the six regions and compared regional differences by calculating Δregional Smax (mean Smax of the highest Smax region – mean Smax of the lowest Smax region) according to concentration or final rhythm (Fig. 5, Supplementary Table 1, Supplementary Fig. 3).Correlating the effects of virtual and clinical AMDWe observed whether the virtual AF persisted for 32 s, including the pacing time. Termination was defined as the absence of an activation signal at the final observation time (32 s). Conversion to AT was defined as the change from AF to regular tachycardia during the observation period (Fig. 1f). Based on the virtual AMD test results, individuals were categorized into the Effective group if AF termination occurred at least once at therapeutic concentrations; otherwise, they were classified into the Ineffective group (Supplementary Fig. 1). Subsequently, we conducted a comparative analysis of the virtual Effective and Ineffective groups to assess the clinical recurrence of AF or AT from the actual date of AMD prescription until 1 year later. A model was created to predict the probability of maintaining sinus rhythm 1 year after AMD treatment.Statistical analysisContinuous variables without a normal distribution are presented as medians and interquartile ranges, and variables with a normal distribution are presented as means ± standard deviations. Continuous variables without a normal distribution were analyzed using the Mann-Whitney U test for two-group comparisons and the Kruskal-Wallis test for three or more group comparisons. Continuous variables with a normal distribution were examined using a t-test for two-group comparisons and analysis of variance to compare three groups.We conducted paired t-tests or Wilcoxon rank-sum tests to assess changes in continuous variables between the virtual AMD test groups. Additionally, we used the Cochran-Armitage test for trend analysis and a linear regression analysis for examining differences in responses based on AMD concentration. A Kaplan–Meier analysis with a log-rank test was used to analyze AF recurrence-free survival over time and compare recurrence rates among the groups. To identify predictors associated with clinical AF recurrence after 1 year of AAD use, a multivariate Cox regression analysis was performed. To assess the predictive ability of maintaining sinus rhythm 1 year after the virtual AMD test, we calculated the area under the receiver operating characteristic curve. Statistical significance was set at a two-sided P-value < 0.05. All statistical analyses were performed using R version 4.2.3 (R Foundation for Statistical Computing, Boston, Massachusetts, United States).
Clinical usefulness of digital twin guided virtual amiodarone test in patients with atrial fibrillation ablation
