Stochastic virtual heart model predictions

Replication of the study "Stochastic virtual heart model predictions" Our study published in "Matters Arising | Nature Cardiovascular Research" (DOI: https://doi.org/10.1038/s44161-025-00641-1) highlights the presence of severe numerical artefacts in the study by Sung et al [1] which uses numerical settings that are incompatible with the used simulation tool. The outcome of the study does not only depend on the used computational model parameters, it strongly and highly sensitively depends on the used spatial resolution, the segementation of the domain in tissue types (healthy versus borderzone) and the fiber orientation relative to the scarred area. Altering any of these factors can toggle the outcome of a simulation, showing the induction of a ventricular tachycardia (VT) circuit that is maintained, or the exact opposite. By altering all three factors in a systematic manner we show that the outcome of the study by Sung [1] is essentially random. In here, we provide scripts used for producing the results of the study, facilitating the replication of essential figures and tables in the main text and in the supplement. Prerequisites for replicating the study Owing to numerous omissions in the methodological description by Sung et al [1], a direct replication of their study or a thorough scrutiny of results is effectively impeded by two major factors: 1) The exact version of the simulator used is ambiguous. The authors claim that the study was performed with the simulator CARP downloaded from our web page at https://carp.medunigraz.at. However, this is incorrect as the simulator CARP was never available for public download from there, nor from any other website or repository. The link served only as a redirect to the web page https://carpentry.medunigraz.at from where a specific version of our proprietaty simulator CARPentry was made available. This software was available in binary form, with a license that expired 2 years prior to the publication by Sung [1]. Inspection of some input files that have been made available within the course of our scientific exchange further corroborate that the simulator specification is wrong. The input files clearly point towards a proprietary version of the simulator CARP implemented by Plank and Vigmond [2], which was never made publicly available. As this software was subject to hundreds to thousands of changes it cannot be determined with certainty which exact version has been used, and from which repository the software has been downloaded. As the numerical underpinnings are, in principle, the same as implemented in the publicly available software openCARP we provide scripts for replicating our findings using the openCARP simulator [3]. 2) None of the essential input data have been made available, effectively precluding the scrutiny of the study. These data include the meshes, the specification of the pacing sites, the pacing protocol used for the induction, or any output datasets that could be used to validate a simulation run. As such, our study uses the descriptions given in the paper, that is, the use of an average resolution of 390µm and the specified pacing protocol, and investigates all issues on geometrically simplified setups where all numerical issues can be more easily highlighted, and do not go unnoticed as in a complex anatomical models where these issues were overlooked by the authors. Software and data for replication include then: openCARP in the release verion v16 as downloaded from here: carputils in the release matching the openCARP simulator forcepss as included in openCARPutils for running the study protocol meshtool for measuring mesh element statistics bash is used for generating command lines and executing the experiments Instructions for performing all numerical experiments For the sake of reproducibility an overall bash script, runme.sh, is provided to produce all the commands required to execute all experiments. Running ./runme.sh --dry prints all individual commands in the required order. Dependence of conduction velocity upon spatial resolution As shown in Fig. 2b) in the main text, conduction velocity depends to a significant extent on spatial resolution. This is particularly problematic for very slow conduction velocities where artificial numerical conduction block is observed. This dependency is shown for the numerical settings and parameters used in the studies due to Arevalo et al [4] and Sung et al [1], referred to as VARP and inFAT settings. The measurement of intrinsic conduction velocitioes along the eigenaxes of the tissue is performed as follows: ./CV_measure.sh -f varp_infat.json Spatial extent of the depolarization wavefront and spatial resolution Numerical conduction block is best understood by measuring the spatial extent of the wavefront. A sound spatial discretization of the reaction-diffusion monodomain equation requires the spatial resolution to be smaller than the spatial width of the depolarization wavefront. The spatial extent of the depolarization wavefront as produced by the VARP and inFAT settings is illustrated in Fig. X in the supplelement. These measurements can be reproduced as follows: ./propagated_AP_inFAT.sh The script generates elongated 3D strand models of high spatial resolution, in the range of 25-50µm, initiates planar wavefront propagation and extracts the spatial distribution of the transmembrane voltage along the strand, Vm(x), at the instant when the wavefront arrives at the center of the strand. The spatial extent of the depolarization wavefront is measured then as the distance between the locations of 20% and 80% action potential amplitude. VARP experiment The scripts for running the entire VARP experiment are provided to replicate our findings. To re-run the entire study, the following series of shell script commands can be used: auto-varp.sh The final simulation results for all tissue discretizations (Nx=6 setups using 50, 100, 200, 300, 350, 400 and 500 µm), fiber orientations (Nf=3 with fibers running at angles 0, 45 and 90 degrees) and two borderzone width (Nbz=2 with 1mm and 2mm) are provided in the form of movies for the sake of saving space (the entire data of the study amount up to about 1 TByte of data) Mesh element statistics In the course of our scientific exchange two patient datasets, P_JHU and P_UPH, have been made available here: GitLab JHU Trayanova For each patient three variants are available, one mesh derived from CT, one from MRI and a fused mesh CTMRI where scar, border zone and fat labels have been registered onto the CT anatomy. In the response (DOI: https://doi.org/10.1038/s44161-025-00642-0) the authors argued that numerical problems have been avoided by using adaptive meshing, where a much finer resolution was used in critical areas of slow conduction. We used meshtool to measure mesh element discretization statistics to verify this statement. Our measurements show that this statement is entirely incorrect. All regions in all six meshes have been meshed around an average resolution of 390µm, there is no difference between normal and slow conducting regions. Our measurements are replicated by running the provided shell script: bash measure_edge_stats.sh A table on the measured average edge lengths in all meshes and all domains is provided below: Table edge length measurements: Model P_JHU_CTMRI Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 393.3 72.4 45.7 675.3 intermediate 2 392.0 72.3 45.7 854.0 scar 3 388.7 74.1 88.8 730.2 inFAT 4 392.8 72.5 45.7 764.1 Model P_JHU_CT Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 393.5 72.3 45.7 691.2 intermediate 2 390.5 73.9 45.7 854.0 scar 3 390.9 73.1 45.7 748.9 Model P_JHU_MRI Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 392.5 73.6 45.2 854.6 intermediate 2 389.5 75.6 45.2 839.0 scar 3 392.8 73.4 45.2 810.7 Model P_UPH_CT Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 393.4 72.5 56.0 699.6 intermediate 2 390.7 73.9 56.0 745.8 inFAT 3 391.1 73.7 56.0 676.0 Model P_UPH_MRI Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 393.9 72.4 56.5 749.2 intermediate 2 391.7 73.0 56.5 683.9 scar 3 392.8 72.5 56.5 667.7 Model P_UPH_CTMRI Tissue type Label Mean Stddev Min Max [µm] [µm] [µm] [µm] normal 1 393.5 72.4 61.2 699.6 intermediate 2 390.8 73.2 61.2 745.8 scar 3 389.4 73.6 66.5 697.2 inFAT 4 391.9 73.2 56.0 682.4 References [1] Sung, E. et al. Fat infiltration in the infarcted heart as a paradigm for ventricular arrhythmias. Nature Cardiovascular Research 1, 933–945 (2022). [2] Vigmond, E. J., Weber dos Santos, R., Prassl, A. J., Deo, M. & Plank, G. Solvers for the cardiac bidomain equations. Progress in Biophysics and Molecular Biology 96, (2008). [3] Plank, G. et al. The openCARP simulation environment for cardiac electrophysiology. Computer methods and programs in biomedicine 208, 106223 (2021). [4] Arevalo, H. J. et al. Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nature communications 7, 11437 (2016). DOI:10.1038/ncomms11437