As in the studies by Martinez-Navarro et al. and Riebel et al. (17, 20), a biventricular model of human electrophysiology embedded in a torso was created through the integration of an anatomical model from clinical magnetic resonance images of a human subject, and extensive multiscale experimental and clinical data from ionic to body surface potentials. Myocardial fibre orientation was represented using a rule-based method (21) reproducing the experimental findings of Streeter et al. (22) (Figure 1A, left). Under healthy conditions, longitudinal conduction velocity along the fibre direction was 65 and 50 cm/s in the transversal and normal directions. The conductivities in the fast conductive endocardial layer were four times higher than in the rest of the ventricles to obtain realistic QRS morphologies and durations.

The human ventricular cell ToR-ORd model (23), with established credibility for healthy, disease, and drug block conditions, was used to represent membrane kinetics. Transmural heterogeneities were defined through a sensitivity analysis to obtain realistic T waves in the electrocardiogram: 50% endocardium, 30% mid-myocardium, and 20% epicardium cardiomyocytes (23). The experimentally reported apicobasal gradient (24) in the slow delayed rectifier potassium current (I_Ks) was also introduced leading to gradients in action potential (AP) duration of 11 ms.

Electrophysiological effects of acute ischaemia, caused by hyperkalaemia, hypoxia, and acidosis, were incorporated into the ischaemic tissue to model the heterogeneous changes in refractoriness and conduction velocity in human ischaemic tissue (17, 25, 26). The formulation of the ATP-dependant potassium current I_K(ATP) (27) was incorporated to simulate hypoxic conditions. Electrophysiological alterations in the ischaemic zone were included as in Martinez-Navarro et al. and Dutta et al. (17, 25), considering the ischaemic core zone (ICZ) and the lateral border zone (BZ), which represents an electrophysiological gradient between the ICZ and the remote myocardium (Figure 1A, right). In addition, the subendocardial border zone described by Wilensky et al. (28) was included to consider oxygen diffusion from blood in the ventricular cavities.

As illustrated in Figure 1B, our simulations included variability in the location of regional ischaemia—left anterior descending artery (LAD), left circumflex artery (LCX), and right coronary artery (RCA) occlusions, see rows in Figure 1B—transmural extent (transmural and subendocardial, see columns in Figure 1B), and ischaemic severity in the remote myocardium due to impaired myocardial circulation (grey scale: healthy tissue, mild ischaemia and moderate ischaemia, described in Figure 1C). ICZ severity was the same in all scenarios (moderate ischaemia), corresponding to the first minutes post-occlusion, as seen in previous studies (17, 19, 26).

To evaluate the impact of ischaemia in the remote myocardium in fibrillatory dynamics, each of the seven scenarios was simulated under three simulation protocols (Figure 1D). In protocol I, 12 s of electrical activity are simulated with healthy remote myocardium. In protocol II, the remote myocardium transitions from healthy to mildly ischaemic after 6 s. Finally, in protocol III, the remote myocardium is mildly ischaemic for 12 s and becomes moderately ischaemic for 6 additional seconds (representing conditions of global myocardial ischaemia).

Sinus rhythm was simulated to obtain realistic ECGs by implementing endocardial stimulation in root nodes (four in the left ventricle, three in the right ventricle) and a fast endocardial layer representing the Purkinje network (29, 30). VF was induced by a stimulation protocol with progressively higher pacing frequency. The state variables of the cellular models in the biventricular simulations were loaded from single-cell simulations with a cycle length (CL) of 350 ms, for each cell type and ischaemic severity. Then, two regular beats at CL = 350 ms were simulated, followed by a sequence of sinus rhythm stimuli with progressively shorter CL (3× CL = 325 ms, 5× CL = 300 ms, 5× CL = 275 ms, 6× CL = 250 ms), taking approximately 6 seconds—last stimulus set when t = 5,850 ms—to complete the pacing protocol (black rectangle in Figure 1D). No further pacing was induced until the end of the simulations to assess VF sustenance.

The ECG signal was computed in four bipolar electrode locations typically considered in external defibrillation protocols: (a) apex-anterior, (b) anterior-lateral, (c) apex-posterior, and (d) anterior-posterior electrodes, as shown in Figure 1E (31). In addition, the ECG was computed at standard clinical precordial lead locations in the 12-lead electrocardiogram (e).

ECG signal during VF was analysed in the frequency domain using the fast Fourier transform (FFT) to obtain the DF and AMSA (10, 19). AMSA was calculated as the sum of the products of each frequency and their amplitudes: AMSA = ∑A_i · F_i, where A_i represents the amplitude at the ith frequency F_i. The calculation parameters reproduced those in the study by Ristagno et al. (10): the signal was resampled at 250 Hz and analysed on a 512-point window (2.05 s), and a Tukey FFT window (α = 0.2) was used to minimise edge effects in the analysed window. An additional ECG marker was calculated, namely, the MS of the signal as in the study by Neurauter and Strohmenger (9). For consistency, DF and MS were measured on the same time windows as AMSA.

Wave dynamics throughout the ventricles in VF were quantified as in Martinez-Navarro et al. (17) by identifying the number of PS, calculated as the lines formed by overlapping two isosurfaces in the ventricles: V_m = −40 mV and dV_m/dt = −2. The number of PS was determined in the midpoint of the time windows used for electrocardiographic markers.

The original formulation of the ToR-ORd human cardiomyocyte model exhibits a relatively shallow action potential duration (APD) restitution curve and limited capability to trigger alternans, both important for fibrillation (32). Therefore, to optimise the model for VF investigations, we created a population of 1,000 variations of the ToR-ORd model using Latin Hypercube Sampling by varying all current conductances and the inactivation dynamics in currents critical for restitution dynamics (I_Na, I_NaL, I_CaL, I_Kr). Firstly, only those models (N = 475) reproducing experimental data in humans describing AP morphology and intracellular calcium measurements (33) were considered (Supplementary Figure S1). A second calibration step was performed to retain only the models that reproduced the AP duration ranges reported in Sutton's work (26) under ischaemia conditions (phase 1A) and a pacing frequency of 2 Hz (Supplementary Figure S2), resulting in 275 variations of the ToR-ORd model. The S1S2 restitution protocol was applied to these models, selecting the 17 variations with the steepest restitution curve. Finally, the dynamic restitution protocol was applied to these models (and baseline ToR-ORd) to investigate their capability to yield alternans. For this study, we chose the ToR-ORd version with the steepest AP restitution and greatest APD alternans amplitude (Supplementary Figures S3A,B, respectively) assuming this model would be the most prone to trigger spiral breakup and alternans. Only models presenting no repolarisation abnormalities at short CLs were considered valid (Supplementary Figure S4). The features of the selected model variant, ToR-ORd_VF, are achieved by a combination of factors, including slower inactivation dynamics in the I_Na (50% decrease of τ_h and τ_j) and I_NaL (34% decrease of τ_h,L) currents, as well as increased conductance of the same currents (39% increase for both). These changes altered tissue propagation, leading to a 5% reduction in the myocardial conduction velocity. A complete characterisation of the model is provided in the Supplementary Methods (Supplementary Figure S5). In addition, tissue conductivities were reduced to yield conduction velocities ∼20% slower than in the healthy myocardium to promote arrhythmia (34, 35), accounting for a 25% conduction velocity reduction compared to the baseline model at tissue level.

All simulations were conducted on the open-source solver MonoAlg3D (36) using graphics processing units (GPUs). For verification, a set of simulations computed in MonoAlg3D was verified against those obtained from other cardiac electrophysiology solvers, such as Chaste (30). All simulations produced identical ECG morphologies and those run on MonoAlg3D incurred lower computational costs and hardware requirements. Simulations were run in the supercomputers Piz Daint (Swiss National Supercomputing Centre), Cirrus (Edinburgh Parallel Computing Centre, UK), and Polaris (Argonne Leadership Computing Facility, USA). About 10 s of simulated human cardiac electrophysiology took ∼18 h, incurring a lower cost than simulations run on solvers that scale only on a multi-CPU paradigm (∼7× speedup).

Model credibility is supported by validation of the original human ventricular cell ToR-ORd model with extensive human data, as described by Tomek et al. (23) for healthy, diseased, and drug block conditions. In the new ToR-ORd_VF formulation, AP morphology and calcium dynamics both in healthy and ischaemia conditions reproduced experimental data in human (26, 33) while exhibiting alternans and a steeper AP restitution. Furthermore, simulations in MonoAlg3D with ToR-ORd_VF showed agreement with experiments in human tissue (37), and ECGs obtained with biventricular simulations were evaluated by comparing the simulated ECGs with clinical records as shown in the next section.

Figure 2 highlights the agreement of simulated and clinical ECGs, supporting the credibility of simulations in healthy (Figure 2A) and ischaemic conditions (Figure 2B). Simulated ECGs using the baseline ToR-ORd and the optimised ToR-ORd_VF cell models show similar patterns to representative clinical ECGs from healthy patients obtained from the study by Mincholé et al. (30) and the UK Biobank (38). T waves in healthy patients show variable amplitudes and morphologies but are typically positive in V3 (0.3–0.6 mV), as shown in Figure 2A (right panel). Figure 2B shows how simulations using the baseline ToR-ORd model can reproduce the distinctive ECG signature of acute ischaemia (39, 40), such as hyperacute T waves and ST elevation ≥0.2 mV (left panel), as in our previous work (17). The simulated ECGs are compared with ECGs from ischaemic patients accessed via the PhysioNet repository (https://physionet.org/), specifically the STAFF database, collected with Institutional Review Board approval (40). Interestingly, the inclusion of the ToR-ORd_VF model characterised by a steep APD restitution (see Section 2.5 for full model description), allowed us to reproduce more complex ischaemia-induced ECG abnormalities in our simulations. Specifically, we were able to reproduce the tombstone pattern, characterised by an elevated ST segment merged with the T wave, larger than the preceding R wave, which tends to be short and of small amplitude (41) (Figure 2B, right panel). This pattern has been observed in certain anterior wall ischaemia/infarction patients (42) and is linked with higher mortality (43).

To further support model credibility at the cellular level, Figure 2C illustrates that both the baseline ToR-ORd and ToR-ORd_VF models (dashed and solid lines, respectively) yield APD₉₀ within experimental ranges in control (grey region) and ischaemic conditions (purple region) (26). Simulations with ToR-ORd_VF produce a steeper and realistic (44, 45) restitution curve (maximal steepness of 0.6 when DI > 280 ms, 3.7 when DI < 280 ms) and increased capability to reproduce alternans while showing no repolarisation abnormalities at physiological pacing frequencies. More information is provided in the Supplementary Methods section.

Simulated ECGs in VF simulations were compared with clinical recordings from the PhysioNet repository (https://physionet.org/), and specifically the Creighton University Ventricular Tachyarrhythmia Database (46, 47), collected with Institutional Review Board approval. Figure 3A shows a representative clinical ECG recording from a patient during a VF episode: after sinus rhythm is replaced by an arrhythmia, DF accelerates progressively from 5 to 8 Hz; 24 seconds after arrhythmia onset, signal amplitude decreases from 2 to 1 mV. These changes, happening within a minute from VF onset, suggest a continuously changing substrate caused by secondary acute myocardial ischaemia (48). Figure 3B presents simulated ECGs during anteroseptal ischaemia leading to VF. We simulated the first few minutes of cardiac activity after arrhythmia onset by increasing the severity of the ischaemic substrate, represented as grey horizontal bars. As a result, we observed similar changes in DF and signal amplitude to those in the clinical record.

For further comparison with clinical and experimental data, we quantified VF dynamics and electrocardiographic markers in five time windows across the simulations for ventricular tachycardia/fibrillation (VT/VF) onset (Figure 3B, window 1, t = 3,905–5,950 ms), VF persistence (Figure 3B, windows 2 and 3, t = 5,950–7,995 ms and t = 9,955–12,000 ms, respectively), and transition to global ischaemic conditions (Figure 3B, windows 4 and 5, t = 12,000–14,045 ms and t = 15,955–18,000 ms, respectively).

Figure 3C further highlights the agreement of simulated dynamics with reported ECG and activation patterns extracted for the classic VF stages defined by Wiggers (48). Wavefronts are shown as black bands over a grey background and PS as magenta circles. •

Sinus rhythm. Simulation of normal myocardial activation at an increasingly fast pacing frequency (CL = 350–250 ms) to induce VF. The ECG presents ST elevation, typically associated with acute myocardial ischaemia.

We quantified arrhythmia inducibility for the different regional ischaemia configurations for healthy versus ischaemic remote myocardium. As shown in Figure 4, three out of seven cases featuring a healthy remote myocardium (left column, simulation protocol I) led to transient VT/VF episodes, as they self-terminated within the 12 seconds simulated. Out of these three VT/VF cases, two of them led to persistent VF if ischaemia was introduced in the remote myocardium after VF induction (centre column, simulation protocol II). Finally, VT/VF occurred in all scenarios with ischaemic remote myocardium, and six out of seven cases presented persistent VF throughout the whole simulation (right column, simulation protocol III).

Figure 5 illustrates the manifestation of VF dynamics in the ECG through three ECG-extracted metrics (AMSA, MS, DF) for five electrode configurations compared to PS. The electrode configurations considered are commonly used in resuscitation protocols: apex-anterior (Figure 5A), normalised apex-anterior for 1 mV amplitude (Figure 5A2), anterior-lateral (Figure 5B), apex-posterior (Figure 5C), and anterior-posterior (Figure 5D) locations.

Each panel includes a graphical representation of the electrode placement, and plots of ECG markers—AMSA, MS, and DF in the left, centre, and right columns, respectively—vs. the number of PS. Circles in the scatter plots represent a measurement from the simulations. The colour of the circles codes the degree of ischaemia in the remote myocardium: white (regional ischaemia + healthy remote myocardium), grey (regional ischaemia + mildly ischaemic remote myocardium), and black (global moderate ischaemia). Data points are visually clustered in black rectangles to represent the value ranges measured for each ischaemic configuration. Measurements obtained from the VT/VF onset stage (circles with orange border) were excluded in the clustering, due to the high disparity in the number of PS and ECG marker values shown during spiral breakup. The red dashed line in AMSA and MS plots shows the minimal threshold values linked with positive patient outcomes in defibrillation (10, 49) to provide the context within reference clinical values to our simulation results.

The number of PS was consistent among panels in Figure 5 as the cases considered for each electrode configuration were the same. VF in cases with a healthy remote myocardium (white points) was uncommon, with only three data points excluding VT/VF onset measurements (1, 4, and 15 PS, white circles). The number of PS stabilises in the range of 11–22 in cases with a mildly ischaemic remote myocardium (grey circles) and 3–12 in global moderate ischaemia (black circles), in line with experimental human measurements (6).

AMSA, MS, and DF showed differences in value ranges depending on both ischaemic conditions and the ECG electrode placement. Apex-anterior electrodes (Figure 5A) were considered representative of the different electrode configurations. AMSA values (Figure 5A, left scatter plot) varied according to the severity of ischaemia in the remote myocardium (healthy: 23.62–24.45 mV Hz; mild ischaemia: 10.58–21.47 mV Hz; moderate global ischaemia: 4.82–11.12 mV Hz), allowing the stratification of VF cases according to ischaemic conditions. Hence, the overlap in AMSA value ranges is minor, represented graphically by the black rectangles.

MS values also decreased as the severity of ischaemia increased, but there was substantial overlap in the value ranges observed across different ischaemic conditions (Figure 5A, centre scatter plot; healthy: 4.13–9.25 mV/s; mild ischaemia: 3.54–5.71 mV/s; moderate global ischaemia: 1.82–4.44 mV/s). Regardless of ischaemic severity, stratification remained challenging around the threshold of MS = 4.2 mV/s. The red dashed line represents the clinical reference value from the study by Shanmugasundaram et al. (49) (2.6 mV/s). Finally, DF values (Figure 5A, right scatter plot) also varied according to the severity of ischaemia in the remote myocardium. In the healthy case, DF oscillates from normal rhythm to tachycardia values (1.5–5.5 Hz). A wide range of DF values was observed in mild ischaemia (3.0–9.0 Hz) but global moderate ischaemia stabilised DF to 7.5–8.0 Hz (except for cases reporting the second harmonic of that frequency, ∼15 Hz).

Subsequently, AMSA was characterised as the most suitable ECG marker for stratifying VF according to the underlying ischaemic substrate. The different electrode locations (Figures 5A–D) showed variability in their discriminatory value. Simulations using apex-anterior (Figure 5A) and apex-posterior (Figure 5C) electrode configurations characterised global myocardial ischaemia cases with AMSA values below the reference value linked to poor resuscitation prognostic, 15.5 mV Hz as in the study by Ristagno et al. (10) (red dashed line). However, the stratification of VF based on the ischaemic substrate in other electrode locations would require different threshold values for AMSA.

While the range of DF values was consistent across electrode locations (1.5–15 Hz), the range of AMSA and MS values changed based on the electrode configuration used. We hypothesised that this change could be due to the variability in the ECG signal amplitude. To confirm this, we compared ECG markers extracted from the same electrode configuration, with and without normalisation. Figure 5A shows ECG markers based on the original simulated ECG signal, and Figure 5A2 presents measurements from the normalised simulated ECG signal (amplitude of 1 mV during VF), which required scaling up the original signal by 23%. DF values remained unaltered post-normalisation (right plot in Figures 5A,B), whereas AMSA and MS measurements increased proportionally to the signal amplitude (left and centre plots, Figure 5A,A2). Even though AMSA can stratify VF cases according to the underlying ischaemic substrate, assessment based uniquely on clinical reference values may not suffice, as the ECG signal amplitude shifts the range of AMSA values adopted.

As shown in Figure 5, ECG biomarkers and number of PS vary with the progression of ischaemia in the remote myocardium, reflecting changes in VF complexity. To provide further mechanistic insights, Figure 6 illustrates this progression through different VF stages: A.

VT/VF onset with healthy remote myocardium (Figure 6A), corresponding to re-entrant patterns with a low range of PS. A single re-entry (1 PS) shown in Figure 6A was sustained through the ischaemic BZ in the septum, forming a re-entrant pathway emerging periodically at the epicardium (5,135–5,400 ms). This activation pattern is reflected in the ECG by a tachycardia-like signature featuring periodic and organised whole-ventricular depolarisations, leading to elevated AMSA (>15.5 mV Hz) and MS values (>2.6 mV/s), linked with a positive patient resuscitation outcome.

Supplementary Video S1 shows a representative example of the VT/VF progression described above.

As shown in Figure 6A and reported extensively in the literature (17, 50), the ischaemic BZ in regional ischaemia provides a powerful substrate for re-entrant wavefronts. In our simulations, regional ischaemia induced the necessary electrophysiological heterogeneities for re-entries leading to tachycardia-like ECG signatures in cases with healthy remote myocardium (see Figure 6A). However, spiral breakup and posterior progression of VF depended on the presence of secondary ischaemia in the remote myocardium. Supplementary Figure S6A shows a VF simulation under global myocardial ischaemia conditions. The membrane potential maps in Panel A represent the fibrillatory wave dynamics (t = 11,760–12,160 ms) at the transition from global mild to moderate ischaemia (first and second row, respectively). The highly disorganised activation pattern based on multiple unstable wavelets (13 PS) quickly evolved into a lower number of very stable rotors (9 PS), and the ECG signature changed accordingly (Supplementary Figure S6B, AMSA dropped from 18.60 to 10.79 mV-Hz). The red rectangle represents the time window visualised in Panel A, and the increase in ischaemic severity is annotated by the blue vertical line. These results suggest that ischaemic severity determined alterations in VF wave dynamics and ECG signature, while heterogeneities caused by regional ischaemia were responsible for the occurrence of re-entrant pathways, potentially leading to arrhythmias.