{"paper_id":"d17c7d1a-eb7c-431a-95f9-d37793671cae","body_text":"MonoAlg3D: Enabling Cardiac Electrophysiology\nDigital Twins with an Efficient Open Source Scalable\nSolver on GPU Clusters\nLucas Arantes Berg1,+,*, Rafael Sachetto Oliveira2,+,*, Julia Camps3, Zhinuo Jenny Wang1,\nRuben Doste1, Alfonso Bueno-Orovio1, Rodrigo Weber dos Santos4, and\nBlanca Rodriguez1\n1Department of Computer Science, University of Oxford, Oxford, United Kingdom\n2Department of Computer Science, Universidade Federal de S˜ao Jo˜ao del-Rei, Brazil\n3Department of Automatic Control, Universitat Polit`ecnica de Catalunya, Spain\n4Graduate Program in Computational Modelling, Universidade Federal de Juiz de Fora, Juiz de Fora, Brazil\n*lucas.arantesberg@cs.ox.ac.uk and sachetto@ufsj.edu.br\n+these authors contributed equally to this work\nABSTRACT\nModelling and simulation are essential in biomedicine, and specifically in computational cardiology. Reliable, efficient and\naccurate solvers are critical. This study presents an open source, GPU-based cardiac electrophysiology solver for scalable\ndigital twin multiscale simulations (MONOALG3D), incorporating conduction system calibration and performance optimization.\nThe solver employs the monodomain equation coupled with the Purkinje network, solved via the finite volume method, featuring\na GPU-based linear solver and concurrent simulation dispatch with MPI. We demonstrate a 10.94× speedup over a CPU-based\nsolution and scalability by running 512 simulations on 128 compute nodes, completing all coarse-mesh simulations in less than\n24 minutes and fine-mesh simulations in 303 minutes. We also demonstrate integration into a cardiac digital twin pipeline for\npersonalisation based on clinical data. The proposed open source solver enhances computational efficiency and physiological\nfidelity, enabling large-scale, high-speed cardiac simulations. This work marks a significant step toward fast and scalable\ncardiac simulations on GPU architectures, with integration in a Digital Twin personalisation pipeline including the conduction\nsystem.\nIntroduction\nComputational modelling and simulation techniques in biomedicine have advanced over the last decades, from enabling\ninvestigations of disease mechanisms to making in silico trials for therapy evaluation possible1–3. Computational cardiology\nis a field that exemplifies a substantial amount of progress 4–7. Credible computer models of the heart are now available from\nsubcellular to whole-organ dynamics, and efficient and accurate solvers have also been made available enabling large simulation\nstudies. These two advances have driven the field forward towards the realisation of the the ‘Digital Twin’ vision in healthcare\n7–10 and in silico trials for therapy evaluation11. For this, cardiac models need to incorporate clinically-relevant features of\ncardiac function and structure, such as the Purkinje conduction system, fibre orientation, anisotropy, cell coupling, and ECG\ncomputation, among others. Moreover, in silico clinical trials and therapy evaluation require consideration of large cohorts of\nvirtual patients. Given ecological and economical limitations in computational resources, simulation software needs to provide\nan accurate and efficient approximation to the mathematical models describing such phenomena. Finally for reproducibility\npurposes, the software and models need to be made open source.\nThe challenge of the high computational costs associated with the resolution of cardiac models has led to the development\nof more efficient numerical schemes and the adoption of parallel computing techniques to reduce simulation times. In addition,\nthe inclusion of the Purkinje system is a fundamental step toward physiological accuracy and correct modelling of ventricular\nactivation in cardiac digital twin models. In experimental and computational studies, it has been shown that this structure can\ninitiate and maintain certain types of arrhythmias due to altered conduction properties under pathological conditions leading\nto ectopic beats and reentrant circuits12.\nEqually important, open source software presents several advantages 13: transparency (as the source code is freely available\nto the public, allowing end-users to validate and verify its functionalities), as well as flexibility and modularity (users can\nadapt and customise the software to their particular needs by adding novel features or implementing new modules that can be\nshared in a collaborative development environment). Along these lines, the cardiac modelling and simulation community has\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\ncontributed several solvers to address these challenges14–16. OpenCARP17 and li f ex-ep18 constitute instances of more recent\nopen source cardiac solvers based on central processing units (CPUs). However, despite providing substantial functionality,\nfeatures such as the simulation of the Purkinje system remain to be addressed.\nThis study presents key improvements to MONOALG3D, an open source, GPU-based solver for cardiac electrophysiology\nsimulations, based on 19. These enhancements significantly improved computational performance and broaden its applicability.\nThe widespread availability of Graphical Processing Units (GPUs) enables substantial speed-ups, which are essential for\nlarge-scale in silico trials 20–24. In addition, GPU-based solvers offer a cost-effective alternative by reducing reliance on\nextensive CPU resources in high-performance computing (HPC) environments. With the increasing adoption of GPU clusters,\ndriven in part by advancements in artificial intelligence, these improvements further support efficient large-scale cardiac\nsimulations.\nMONOALG3D uses the finite volume method (FVM) to simulate the monodomain model on GPU and/or CPU hardware,\nand OpenMP and NVIDIA CUDA to respectively parallelise CPU/GPU computations. Advancements include, firstly, a fully\nintegrated Purkinje network model that allows for retrograde propagation, enabling investigations on its potential role in\npromoting and sustaining complex arrhythmias. Secondly, a novel solver for the diffusion linear system is implemented directly\non GPUs to reduce the computational time of this component. Thirdly, a new output format is included to optimise disk space\nand access. Finally, the framework is expanded to support parallel dispatching across exascale HPC infrastructures, thereby\nfacilitating large-scale, high-throughput simulations studies essential for digital twin and personalised medicine applications.\nPerformance improvements are evaluated on different hybrid CPU/GPU combinations for the device and on two human-\nbased cell models of increasing complexity, as well as compared to space adaptivity features presented in previous work19.\nFinally, to demonstrate the full capabilities of the proposed solver and verify its scalability on GPU clusters under more realistic\nscenarios, our last experiment presents a cardiac digital twin application considering a biventricular simulation with the Purkinje\nsystem and ECG recordings.\nMethods\nMonodomain model\nThe monodomain model is commonly used to describe electrical propagation due to its lower computational cost compared\nto the bidomain model22, 25. In the next equations we present the mathematical models for the myocardium and the Purkinje\nsystem, along with their coupling. We use subscripts P for the Purkinje domain, M for the myocardium domain, and d ∈ {P,M}\nfor the full domain.\nβ\n\u0012\nCm\n∂Vd\n∂t + Iiond (Vd, ⃗ηd)\n\u0013\n= ∇ · (σd∇Vd) +βIstimd (1)\nin Ωd × (0,T ),\n∂ ⃗ηd\n∂t = fd(Vd, ⃗ηd) (2)\nin Ωd × (0,T ),\n(σM∇VM) · ⃗nM = (σP∇VP) · ⃗nP (3)\non ∂ Ωd × (0,T ),\nVd(Xd,0) = Vd,0(Xd), ηd(Xd,0) = ηd,0(Xd) (4)\nin Ωd,\nwhere Vd is the transmembrane potential of either domain, Iiond the total ionic current associated to the cellular model that\ndepends on state variables ⃗ηd, fd the non-linear system of equations encapsulating the dynamics of the state variables, β the\nsurface-to-volume ratio, Cm the membrane capacitance, σd the domain conductivity tensor, and Istimd an external stimulus. The\nmodel is further closed with appropriate Neumann boundary conditions to ensure flux continuity between the myocardium\nand Purkinje domains as given by equation (3), where ⃗nd is the normal vector of the myocardium or Purkinje domain surfaces,\n∂ Ωd. For myocardial surface nodes not coupled to the Purkinje system, equation (3) simply reduces to a standard non-flux\nboundary condition. Initial conditions are provided by equation (4). For the Purkinje domain we consider the one-dimensional\nform of equation (1), while for the myocardium domain its three-dimensional formulation.\nCardiac tissue is known to be comprised of strongly coupled fibres with anisotropic conduction properties. Such fibres\nare defined for each myocardial element by three orthonormal vectors (⃗f ,⃗s,⃗n), where ⃗f lies on the local fibre or longitudinal\ndirection, ⃗s on the sheet or transversal direction, and ⃗n on the normal direction to the fibre. Moreover, associated with each of\nthese vectors, there exist conductivity values σ f , σt, and σn, jointly defining the myocardial conductivity tensor as:\nσM = (⃗f ⊗ ⃗f )σ f + (⃗s ⊗⃗s)σt + (⃗n ⊗⃗n)σn. (5)\n2/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nFinite volume method applied to the monodomain model\nA common technique to efficiently solve the monodomain model is to divide its reaction and diffusion parts using the Godunov\noperator splitting26. Applied to equations (1)–(2), this leads to the solution of two separate problems: a non-linear system of\nordinary differential equations (ODEs)\n∂Vd\n∂t = 1\nCm\n\u0002\n−Iiond (Vd, ⃗ηd) +Istimd\n\u0003\n, (6)\n∂ ⃗ηd\n∂t = fd(Vd, ⃗ηd), (7)\nand a parabolic linear partial differential equation (PDE)\nβCm\n∂Vd\n∂t = ∇ · (σd∇Vd). (8)\nWithin the different numerical techniques available, the FVM offers a robust approach for solving the monodomain model\ndue to its foundation on conservative principles and applicability to diverse geometries 27. This technique discretises the\ncomputational domain into control volumes. Each control volume is associated with a variable of interest, and the governing\nequations are applied to ensure the conservation of this variable across the control volume faces.\nCell Model\nFor the solution of the cellular electrophysiology model described by Eqs. (6,7) M ONOALG3D offers different techniques for\nintegration. It supports both the explicit Euler method as well as Rush-Larsen or other methods based on the generalization of\nmatrix exponential, such as the Uniformization approach28.\nMyocardium modelling\nTo spatially discretise the diffusion term in equation (8), we consider the relation:\nJd = −σd∇Vd, (9)\nwhere Jd (µA/cm2) represents the density of intracellular current flow.\nApplying the divergence theorem and using equation (8), it yields:\nβCm\nZ\nΩd\n∂Vd\n∂t dv = −\nZ\n∂ Ωd\nJd · ⃗nd ds, (10)\nwhere ⃗nd represents the normal vector to the domain surface. This equation is the basic term for deriving the linear system\nof equations associated with the linear PDE.\nWe now particularise the FVM equations for the myocardium. For simplicity, let us consider a tridimensional uniform mesh,\nconsisting of hexahedra with a space discretisation hM. Located at the centre of each myocardial volume (i, j,k) is a node\nwith the transmembrane potential VM as the associated variable of interest. Assuming that the volumetric membrane current\nrepresents an averaged value in each hexahedron, and using (10), we then have:\n\u0012\nβCm\n∂VM\n∂t\n\u0013\f\f\f\f\f\n(i, j,k)\n=\n−\nR\n∂ ΩM JM · ⃗nM ds\nh3\nM\n. (11)\nTo support spatially varying fibre orientation and anisotropy of the myocardial conductivity tensor, the surface integral\ncalculations in equation (11) consider the total sum of flows on the 6 faces of the control volume (each with face area h2\nM) over\na 27-neighbours stencil. This gives:\nh3\nMβCm\n∂VM\n∂t = h2\nM\n6\n∑\nl=1\nJl. (12)\nEach Jl in equation (12) is implicitly calculated by evaluating the spatial derivatives ofVM via second-order finite differences at\ntimestep n + 1, and computing the average conductivity tensor given by(5) at the surfaces of the discretised volume. Altogether,\nthe previous steps lead to a linear system to solve the diffusion equation using the backward Euler method. Additional details\nare provided in Supplementary Material section A.1.\n3/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nFigure 1. Illustration of the three possible configurations for a Purkinje control volume and the Purkinje coupling model. In\npanels (A)-(C), the control volume for which we are calculating the fluxes is depicted in grey. (A) Normal case, where a\nPurkinje control volume is associated with no branch. (B) Branching case, where a Purkinje control volume is linked to Nbi f f\nother Purkinje control volumes. (C) Terminal case, where a Purkinje control volume is coupled to NPMJ myocardium control\nvolumes from the myocardium domain by a fixed resistance RPMJ. Labels P and M indicate the domain where each control\nvolume is located. (D) Simple Purkinje network with single bifurcation, coupling a Purkinje terminal to its five closest\nmyocardium control volumes (coloured in grey). (E) Direction of the flux JPMJ for anterograde (JPMJA) and retrograde (JPMJR)\npropagation, respectively.\nPurkinje modelling\nTo model the Purkinje system, we consider the one-dimensional form of the linear PDE given by equation(8), with time and\nspace discretisations following an equivalent approach to the myocardial case presented above.\nHowever, in the case of a Purkinje control volume, we have to consider three different possible configurations in our\nPurkinje networks (normal, branching, or terminal) as shown in Figure 1. The total flux given by equation (10) is:\nNormal : Jtot = Jxi+1/2 − Jxi−1/2, (13)\nBranching : Jtot =\nNbi f f\n∑\nj=1\nJx j − Jxi−1/2, (14)\nTerminal : Jtot = JPMJ − Jxi−1/2, (15)\nwhere Nbi f f is the number of Purkinje control volumes linked to the bifurcation, and JPMJ is the flux associated at the Purkinje-\nmuscle-junctions. Following a similar approach to the three-dimensional case, the fluxes Jxi+1/2, Jxi−1/2, and Jx j are calculated\nvia finite differences at timestep n + 1 alongside the Purkinje conductivity σP associated to the surface of the discretised\nPurkinje control volume using harmonic means.\nPurkinje–myocardium coupling\nTo model the coupling between Purkinje–myocardium domains, we consider an additional fluxJPMJ. The electrical stimulus\ncoming from the Purkinje system reaches the myocardium at specialised sites called Purkinje-muscle junctions (PMJs), spread-\ning in the endocardium by a distance of approximately 1 mm between each other29. Importantly, PMJs are known to exhibit\na characteristic asymmetric conduction delay due to electrotonic interactions of around 4 − 14 ms on the anterograde direction\n(Purkinje-to-myocardium)30 , and of about 2 − 4 ms when propagation occurs in the retrograde direction (myocardium-to-\nPurkinje)30. This behaviour is characterised as a source-sink mismatch phenomenon since a single Purkinje terminal may need\nto activate a bulk of myocardium tissue31, 32.\nTypically, the PMJ coupling is modelled by a fixed resistance, linking a Purkinje element to several myocardium elements33.\nWe follow this approach by modelling the flux JPMJ (µA/cm2) using a fixed resistance RPMJ and by coupling a single Purkinje\ncontrol volume to its NPMJ closest myocardium control volumes, as shown in Figure 1D. Moreover, the PMJ flux is given as a\nnon-homogeneous Neumann boundary condition by:\n4/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nJPMJ = 1\nh2\nP\nNPMJ\n∑\nk=1\n(VP −VMk )\nRPMJ\n, (16)\nwhere the sign of the flux determines if JPMJ exerts its action in the anterograde or retrograde direction (see Figure 1E).\nNumerical scheme\nFor the iterative solution of the coupled model, we start by solving the reaction terms describing the Purkinje and myocardium\ncellular models, given by the non-linear systems of ODEs in equations (6)-(7). We consider here the forward Euler method for\nsimplicity, albeit MONOALG3D is equipped with more advanced ODE schemes (such as Rush-Larsen and adaptive forward\nEuler). This gives:\nCm\nV n+1/2\nd −V n\nd\n∆t =\n\u0002\n−Iiond (V n\nd , ηn\nd ) +Istimd\n\u0003\n, (17)\nηn+1/2\nd − ηn\nd\n∆t = fd(V n\nd , ηn\nd ). (18)\nPMJ fluxes are computed next based on equation (16), as:\nJn+1/2\nPMJ = 1\nh2\nP\nNPMJ\n∑\nk=1\n\u0010\nV n+1/2\nP −V n+1/2\nMk\n\u0011\nRPMJ\n. (19)\nThe diffusion terms of the Purkinje and myocardium domains involves the solution of the linear system:\nβCm\nV n+1\nd −V n+1/2\nd\n∆t = ∇ · (σd∇V n+1\nd ) +Jn+1/2\nPMJ h2\nP. (20)\nObserve that by computing JPMJ at time n + 1/2, we decouple the Purkinje and Myocardium domains. This enhances the\nsolver’s modularity, allowing different classes to be used for each domain, but at the cost of numerical stability. While a fully\nimplicit solution of the PDE is unconditionally stable, decoupling the two domains results in a conditionally stable scheme,\nwhere RPMJ constrains the maximum time step.\nECG calculations\nAn approximation for the ECG can be computed by assuming that the tissue is immersed in an unbounded volume conductor34.\nThe surface potential can be then calculated using the equation:\nφe = 1\n4πσb\nZ\nΩ\nβIm\n∥⃗r∥ dΩ, (21)\nwhere σb is the bath conductivity, and⃗r is the distance vector between source and field points, the latter essentially the electrode\npositions of the virtual ECG leads. The source term βIm is given by the solution of the diffusive term ∇ · (σ ∇Vm), which is\navailable in every timestep.\nTo efficiently implement this new functionality in MONOALG3D, we implemented the calculations of the equation (21)\nusing OpenMP in CPUs or CUDA on GPUs environments.\nPerformance efficiency strategies\nSolving diffusion on GPUs\nIn previous work19, the linear system linked to the diffusion term in equation (8) was exclusively solved in the CPU using\nan OpenMP version of the conjugate gradient (CG) method.To enable the solution of large linear systems on GPUs, we first\nconverted its sparse matrix representation from the ALG format to a Compressed Sparse Row (CSR) data structure compatible\nwith the cuSparse library. This allows to directly solve the CG on the GPU by using this data structure together with the\nmethods implemented in the cuBLAS library. It is worth noting that the biconjugate gradient (BCG) method is also available\nin this new version.\n5/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nNew output format\nTo minimise disk space usage and improve output performance, we provide novel support for EnSight files as new output\nformat. This format is also compatible with multiple visualisation tools, such as Paraview, allowing most post-processing\nworkflows to be kept unchanged.\nMPI batch\nFinally, for sensitivity analysis and uncertainty quantification studies, MONOALG3D provides a novel feature for the concurrent\ndispatch of multiple simulations using the message passing interface (MPI) standard. Given a baseline simulation and a range\nof parameters, the solver generates automated configuration files for all possible combinations of input parameters. Each\nconfiguration file is then dispatched in parallel using MPI. This allows to upscale more efficiently the number of jobs running\nin HPC environments, enabling for instance to perform hundreds of simultaneous simulations for a given patient using a wide\na range of parameters. Such a feature is an important step towards in silico trials, drug therapy, and risk assessment studies.\nComputational simulations\nTwo sets of experiments were used to evaluate the improvements implemented in MONOALG3D: a benchmark cuboid mesh\nto quantify performance improvements; a human biventricular mesh coupled to a Purkinje network (see Supplementary Figure\nS4A section A.3), as an exemplar of cardiac digital twin application.\nAll our numerical experiments were performed in the Polaris supercomputer provided by the Argonne Leadership Computing\nFacility, a 560 node HPE Apollo 6500 Gen 10+ system. Each computing node is equipped with a 2.8 GHz AMD EPYC Milan\n7543P 32 core CPU with 512 GB of DDR4 RAM and four NVIDIA A100 GPUs.\nBenchmark myocardium cuboid\nA numerical test, adapted from Niederer et al.35, was conducted with minor domain size modifications to: 1) evaluate GPU\nspeedups in solving non-linear ODEs and the parabolic PDE; 2) assess space adaptivity effects on efficiency; and 3) compare\ndisk space usage between EnSight and VTK formats.\nThe test used a 1 × 1 × 1 cm3 myocardium cuboid with transverse anisotropic conduction σ∥ = 1.334 mS/cm, σ⊥ =\n0.176 mS/cm), monodomain parameters β = 1400 cm−1, Cm = 1 µF/cm2, and human-based ventricular models: ten Tusscher\n(12 state variables)36 and ToR-ORd (43 state variables)37. Stimulation was applied in a 0.15 × 0.15 × 0.15 cm region for 2 ms\nat 53 pA/pF.\nThe nonlinear ODEs were solved using the Rush-Larsen scheme 38 with ∆t = 0.01 ms, while the parabolic PDE used\n∆t = 0.02 ms for a total of 1500 ms. Without space adaptivity, uniform discretization at hM = 250 µm led to 64,000 control\nvolumes. Adaptive resolutions ranged from hMmin = 250 µm to hMmax = 500 µm, with refinement/de-refinement bounds at\n10.01/10.00 and adaptation every 10 timesteps.\nThe benchmark was tested across six CPU/GPU configurations:\n• A+OC+PC: Adaptive, ODEs on CPU, PDE on CPU;\n• A+OG+PC: Adaptive, ODEs on GPU, PDE on CPU;\n• OC+PC: Non-adaptive, ODEs on CPU, PDE on CPU;\n• OC+PG: Non-adaptive, ODEs on CPU, PDE on GPU;\n• OG+PC: Non-adaptive, ODEs on GPU, PDE on CPU;\n• OG+PG: Non-adaptive, ODEs on GPU, PDE on GPU.\nMesh geometry and transmembrane potential were saved every 100 timesteps. More details on the benchmark setup are\nprovided in Supplementary Figure S3 section A.2.\nCardiac digital twin with Purkinje network\nTo demonstrate the full capabilities of the proposed GPU cardiac solver, we conducted a simulation study within a biventricular\ncardiac digital twin pipeline incorporating a Purkinje network. The study had three objectives: 1) evaluate solver performance\nin realistic scenarios; 2) calibrate Purkinje coupling parameters RPMJ and NPMJ to physiological anterograde PMJ delays; and\n3) verify solver scalability for concurrent GPU simulations.\nWe used a human biventricular mesh (76-year-old female,87 kg, 107 cm3 volume) reconstructed from MRI39, previously\napplied in clinical ECG personalization 40–42. Supplementary Figure S4 section A.4, presents further anatomical details,\nincluding its Purkinje network coupling (Supplementary Figure S4A), fiber orientation field (Supplementary Figure S4B),\n6/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nsubendocardial Purkinje coupling layers (Supplementary Figure S4C), and IKs scaling factor map for T-wave personalization41\n(Supplementary Figure S4D).\nFor cellular electrophysiology, we used the ToR-ORd human-based ventricular model37 with modifications for T-wave\npersonalization42: 50% IKr scaling, 5× IKs scaling43, and reducing τ jca from 75 to 60 ms. The Purkinje domain was modeled\nwith the human-based PurkinjeTrovatomodel44. Both ODE systems were solved via the Rush-Larsen scheme with∆t = 0.01 ms.\nPDEs used the same discretization step for a total simulation time of 600 ms. The stimulus protocol consisted of a single\npulse applied at the His bundle (Ncells = 25) with 40 pA/pF amplitude and 2 ms duration. For more details about the mesh\nconfiguration refer to Supplementary Material sections A.3 and A.4.\nTo evaluate performance, we tested two myocardium space discretizations: a coarse mesh (hM = 500 µm) with 855,670\ncontrol volumes and a fine mesh (hM = 250 µm) with 6,845,360 volumes. The Purkinje domain used a fixed hP = 250 µm\nwith 7,948 volumes. Solver scalability was assessed across three simulation setups:\n• 1N1S: 1 node, 1 simulation;\n• 1N4S: 1 node, 4 concurrent simulations;\n• 128N512S: 128 nodes, 512 concurrent simulations.\nA large-scale simulation study calibrated RPMJ and NPMJ within a physiological range using MONOALG3D’s MPI batch\nprocessing (see Supplementary Material section A.5). We ran 512 concurrent simulations, varying RPMJ and NPMJ based on\nan initial calibration (see Supplementary Material section A.6). For the coarse mesh, RPMJ spanned [100,1300] kΩ (32 values)\nand NPMJ [15,50] (16 values). For the fine mesh, RPMJ ranged from [500,2300] kΩ. ECG comparison with clinical data was\nperformed using Pearson’s correlation coefficient across all 8 leads (I, II, V1–V6).\nResults and discussion\nMyocardium cuboid benchmark\nAn initial test was conducted using the OC+PC configuration, which uses entirely the CPU without space adaptivity to solve\nboth the ODE and PDE systems, to evaluate the optimum number of OpenMP threads for the selected HPC facility. Five\nsimulations were executed per number of threads, considering the human ventricular cellular ToR-ORd model(Fig. 2A). The\nbest total execution times were found for an optimal number of 8 OpenMP threads, leading to a ≈ 6.66× efficiency speed-up,\nand enabling benchmark execution times around 30 minutes.\nSimilarly, input/output efficiency was optimised by considering5 broadly adopted scientific formats: VTK-text (ASCII),\nVTK-binary, VTK-binary-compressed, EnSight-text (ASCII), and EnSight-binary (Fig. 2B). The results from this analysis\nhighlight substantial savings in output file size when saving each model state variable (transmembrane potential in our case) in\nEnSight-binary format. This resulted in file storage sizes of merely0.19 gigabytes, while VTK-text required around6 gigabytes\nto store the same outputs. Therefore, the EnSight-binary can save approximately 31×, 25×, 3.5× and 3.23× more disk space\nwhen compared to VTK-text, VTK-binary, VTK-binary-compressed, and EnSight-text formats, respectively.\nWe then evaluated the solver’s performance for each of the 6 considered combinations of CPU/GPU architectures (see\nsection Computational simulations/Benchmark myocardium cuboid), using either the ten Tusscher (Fig. 2C) or the ToR-ORd\n(Fig. 2D) cellular models. Based on the analysis above, 8 OpenMP threads were used in all the cases. The efficiency results\npresented in Figure 2C for the ten Tusscher model yielded a maximum simulation time of ≈ 9 min for the OC+PC scenario\n(i.e., solving the entire problem in the CPU, without space adaptivity). Space adaptivity allowed the problem to be solved under\n6 min in the CPU (A+OC+PC scenario), and in around 4 min if the ODE system was solved in the GPU (A+OG+PC scenario).\nThe largest efficiency improvement was however found when the problem was solved entirely in the GPU (OG+PG scenario),\ndecreasing the simulation time under 2 min. Equivalent results are presented in Figure 2D for the ToR-ORd model for all the\nCPU/GPU configurations. Solving the simulation entirely on the CPU without space adaptivity (OC+PC scenario) yielded\nthe most demanding execution time of ≈ 25 min, while a 10.94× efficiency gain and a total simulation time below 3 min were\nattained by exploiting the full GPU implementation (OG+PG scenario).\nThe results above indicate that space adaptivity (scenarios A+OC+PC and A+OG+PG) did not lead to any improvements\nin performance when compared to solving the fully refined mesh entirely on the GPU (scenario OG+PG). This behaviour can\nbe attributed to the computational overhead associated with spatial adaptivity, specifically reassembling the matrix of the PDE\nand updating the grid data structures. In contrast, preloading the fully refined mesh and transferring all the data structures for\nboth the ODEs and PDE to the GPU at the onset of a simulation offers significant advantages in terms of memory usage and\ncomputational performance when a fixed spatial discretisation is used. This approach not only reduces data transfer between\nthe host and the GPU, but also minimises memory allocation operations. Consequently, the computation at each time step is\nmore regular than when space adaptivity is used, leading to an improved overall performance.\n7/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nFigure 2. Results for the myocardium cuboid benchmark. (A) Execution time of the ToR-ORd simulation considering an\nOC+PC solver. (B) Disk usage for storing the mesh geometry and with the myocardial transmembrane potential for the\nToR-ORd simulation using different file formats. (C) Total execution time for each of the 6 scenarios when using theten\nTusscher model. (D) Total execution time for each of the 6 scenarios with the ToR-ORd model.\nIn addition, the joint analysis of the two considered cellular models revealed that solving the ODE system on the GPU was\nresponsible for most of the performance gains. This becomes apparent by comparing scenarios OC+PC and OC+PG. For the\nten Tusscher model, scenario OC+PG is just 1.11× faster than scenario OC+PC, while for ToR-ORd a speedup of 1.04× is\nobtained. However, when the ODE system was solved on the GPU (scenario OG+PC), a 3.26× gain was attained for the ten\nTusscher model when compared to scenario OC+PC, while this improvement was even more pronounced for ToR-ORd, and\naround 6.60× gain. This result is primarily underlain by differences in algebraic complexity between both cellular models: the\nten Tusscher model consists of 12 state variables, compared to 43 in ToR-ORd. Additionally, the ToR-ORd model involves a\ngreater number of algebraic expressions, making it well-suited for GPU-based computations. Detailed execution times for both\nmodels can be found in the Supplementary Tables S1 and S2 section A.2.\nCardiac digital twin with Purkinje network\nBased on the results of our previous benchmark, all cardiac digital twin simulations were executed entirely on the GPU using 8\nOpenMP threads, exploiting one GPU device for both the Purkinje and myocardium domains as this configuration demonstrated\nthe best performance.\nFigure 3 illustrates selected simulation results at fine mesh resolution for the identified optimal set of Purkinje coupling\nparameters from a total of 512 executions. This set was chosen based on similarity between the simulated and clinical ECG\nsignals across all leads, as well as replicating a range of physiological anterograde PMJ delays across all Purkinje terminals\nwhen activating the myocardium. Nevertheless, the existence of other combinations of Purkinje coupling parameters yielding a\nsimilar activation pattern (Figs. 3A–3B) indicates that a range of cardiac digital twins can be approximated by the considered\nparameters, RPMJ and NPMJ. Ventricular activation started ≈ 48 ms after His-bundle pacing (Figs. 3C–3D), with the whole\nbiventricular domain being activated in around 110 ms (end of simulated QRS complex, Fig. 3E), also within the physiological\nrange of 80 − 120 ms for healthy subjects29. Moreover, the inclusion of heterogeneity in IKs from42 generated a reasonable\napproximation for the T-wave (Fig.3E), with an average PCC of 0.81 across all the 8 independent ECG leads (I, II, V1–V6).\nThe effects of anterograde PMJ delays were correctly recovered as shown in Figure 3D, as well as in Figure 3F for a PMJ\nsite on the right ventricle. As it can be seen in the latter, the closest myocardial control volume to the Purkinje terminal could\nnot be activated instantaneously, due to the source-sink mismatch between the Purkinje terminals and myocardial cells31, 32, 45.\nThe myocardial control volume only became entirely depolarised approximately 4.34 ms after the stimulus reached the PMJ\n8/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nFigure 3. Results for the 512 digital twin simulations with the fine mesh. (A) Average anterograde PMJ delay across all\nPurkinje terminals. (B) Average PCC across all leads between the clinical and simulated ECG. Selected Purkinje coupling\nparameters (RPMJ = 1777 kΩ,NPMJ = 38) are highlighted by red squares. (C-D) Selected cardiac digital simulation with an\naverage PCC of 0.81 in ECG reconstruction and average anterograde PMJ delay of 2.28 ± 2.45 ms, at times t = 40 and\nt = 50 ms, respectively. (E) Comparison between clinical and simulated ECGs. (F) Action potential upstrokes for the PMJ site\nhighlighted in panel (D) with anterograde PMJ delay of 4.34 ms, at the terminal Purkinje volume (blue) and its closest coupled\nmyocardium volume (lime).\nTable 1. Execution times (minutes) for the cardiac digital twin test for the 3 submission scenarios. Total : total time; W rite:\nwriting time; ECG: ECG computation time; MODE/MPDE: time to solve the myocardium ODE/PDE system; PODE/PPDE: time\nto solve the Purkinje ODE/PDE system; MPIend: time for the MPI process to finish.\nMesh Scenario Total W rite ECG M ODE MPDE PODE PPDE MPIend\n1N1S 15.15 0.08 0.07 7.78 4.27 0.20 1.01 -\nCoarse 1N4S 22.24 ± 0.21 0.07 ± 0.00 0.11 ± 0.00 10.94 ± 0.11 7.23 ± 0.09 0.23 ± 0.00 1.05 ± 0.01 23.15\n128N512S 22.41 ± 0.32 0.07 ± 0.00 0.11 ± 0.02 11.06 ± 0.17 7.24 ± 0.15 0.23 ± 0.01 1.06 ± 0.01 23.55\n1N1S 221.49 0.52 0.42 97.62 87.70 0.26 1.13 -\nFine 1N4S 210.66 ± 0.32 0.50 ± 0.00 0.42 ± 0.00 93.80 ± 0.18 84.81 ± 0.05 0.39 ± 0.01 1.13 ± 0.00 213.13\n128N512S 211.94 ± 7.33 0.49 ± 0.01 0.42 ± 0.00 93.87 ± 0.50 85.89 ± 5.56 0.37 ± 0.03 1.13 ± 0.01 302.29\nsite, exhibiting an initial brief spike followed by a distinct blunted depolarisation (Fig. 3F). This behaviour, reported in different\nexperimental studies46, 47 and attributed to the electrotonic effects at the junctions, is further explored in the PMJ calibration\nresults presented in Supplementary Material section A.6.\nIn terms of scalability, Table 1 summarises execution times for our considered submission scenarios. For the 1N1S scenario\n(1 computing node, 1 simulation), total execution times for the coarse and fine discretisations were around 15 and 221 min,\nrespectively. Similar execution times without significant performance loss were observed when4 concurrent simulations were\nexecuted in the same compute node (scenario 1N4S), with average total times around 22 and 210 min for the coarse and fine\nmeshes, respectively. A similar behaviour was observed using the MPI batch feature for the 128N512 scenario, with execution\ntimes of approximately 22 and 212 min for coarse and fine resolutions, respectively. Considering the times for the MPI process\nto start and end, these values were around 23 and 302 min, respectively. From the results in Table1, it also transpires that the\nmost demanding component is the solution of the ODE system of the myocardium, contributing around 11 and 94 min for the\ncoarse and fine mesh resolutions, respectively. This is explained due to the larger number of control volumes of this domain\ncompared to the Purkinje one, making this section of the problem more computationally demanding and the overall bottleneck.\nBased on these results, the novel MPI batch feature illustrates the scalability and efficiency of the proposed cardiac solver to\nconduct human-based cardiac digital twin studies under cluster GPU environments, enabling the correct adjustment of sensitive\n9/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nPurkinje coupling parameters to physiological ranges for anterograde PMJ delays. Additional results on the performed cardiac\ndigital twin simulations are presented in Supplementary Material section A.7. Of importance, they provide further evidence\nsupporting the existence of multiple possible cardiac digital twins with a similar ECG (see Supplementary Figures S9 to S11\nin section A.7, for example). These results also indicate that, while sustaining analogous ECGs, different combinations of\nPurkinje coupling parameters can generate distinct distributions of PMJ delays across the Purkinje terminals. Based on that,\ncertain Purkinje terminals exhibit more variability in their associated PMJ delays, which might indicate a more important role\nof such PMJs to the whole ventricular activation. Another relevant finding from these simulations is the existence of an almost\nlinear relation between the Purkinje coupling parameters and the appearance of propagation block, as analysed by Figures 3A,\nS7 (Supplementary Material section A.6) and S9A (Supplementary Material section A.7).\nConclusion\nIn this work, we have presented an open source, high-performance GPU solver for electrophysiology simulations. By\nsystematically evaluating the solver’s performance across various CPU/GPU configurations, we demonstrated its efficiency\nin solving the monodomain model under different scenarios. Specifically, our findings highlight that leveraging GPUs for\nboth the non-linear system of ODEs and the parabolic PDEs significantly accelerates computations, particularly for complex\nhuman-based cellular models. These results contrast with our earlier study emphasizing space adaptivity, as we found that\ndirectly copying all required data structures to the GPU at the start of the simulation yields greater efficiency gains. Furthermore,\nthe solver’s integration into a cardiac digital twin pipeline demonstrated its scalability in a GPU cluster environment. The ability\nto execute 512 simulations concurrently across 128 compute nodes, with execution times comparable to single-node setups,\nshows its robustness for large-scale studies. These simulations were instrumental in calibrating Purkinje coupling parameters to\nachieve physiological anterograde PMJ delays, showcasing the solver’s applicability in patient-specific modelling that consider\nthe cardiac conduction system as an essential component of the model. Our results affirm the solver’s capability to perform\nlarge-scale monodomain simulations efficiently, including detailed modelling of the Purkinje-muscle-junctions, a feature that, to\nthe best of our knowledge was not implemented in any other open source cardiac solver. Furthermore, the experiments presented\nin this work can be seamlessly expanded to a large cohort of virtual patients, which is a relevant step towards applicability\nin in silico clinical trials and therapy evaluation. This work provides a significant advancement in the field, offering an open\nsource scalable tool for researchers to explore complex cardiac electrophysiology scenarios with remarkable efficiency.\nReferences\n1. Jean-Quartier, C., Jeanquartier, F., Jurisica, I. & Holzinger, A. In silico cancer research towards 3R. BMC cancer 18, 1–12\n(2018).\n2. Delp, S. L. et al. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE\ntransactions on biomedical engineering 54, 1940–1950 (2007).\n3. Musuamba, F. T. et al. Scientific and regulatory evaluation of mechanistic in silico drug and disease models in drug\ndevelopment: Building model credibility. CPT: Pharmacometrics & Syst. Pharmacol. 10, 804–825 (2021).\n4. Fenton, F. H., Cherry, E. M., Hastings, H. M. & Evans, S. J. Multiple mechanisms of spiral wave breakup in a model of\ncardiac electrical activity. Chaos 12, 852–892 (2002).\n5. Lopez-Perez, A., Sebastian, R. & Ferrero, J. M. Three-dimensional cardiac computational modelling: methods, features\nand applications. Biomed. engineering online 14, 1–31 (2015).\n6. Dasí, A. et al. In-silico drug trials for precision medicine in atrial fibrillation: From ionic mechanisms to electrocardiogram-\nbased predictions in structurally-healthy human atria. Front. Physiol. 13, 966046 (2022).\n7. Oliveira, R. S. et al. Ectopic beats arise from micro-reentries near infarct regions in simulations of a patient-specific heart\nmodel. Sci. reports 8, 16392 (2018).\n8. Li, L. et al. Towards enabling cardiac digital twins of myocardial infarction using deep computational models for inverse\ninference. IEEE Transactions on Med. Imaging (2024).\n9. Corral-Acero, J. et al. The ‘Digital Twin’ to enable the vision of precision cardiology.Eur. Hear. J. 41, 4556–4564 (2020).\n10. Niederer, S. A., Sacks, M. S., Girolami, M. & Willcox, K. Scaling digital twins from the artisanal to the industrial. Nat.\nComput. Sci. 1, 313–320 (2021).\n11. Viceconti, M. et al. In silico trials: Verification, validation and uncertainty quantification of predictive models used in the\nregulatory evaluation of biomedical products. Methods 185, 120–127 (2021).\n10/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\n12. Haissaguerre, M., Vigmond, E., Stuyvers, B., Hocini, M. & Bernus, O. Ventricular arrhythmias and the His-Purkinje\nsystem. Nat. Rev. Cardiol. 13, 155 (2016).\n13. Kavanagh, P. Open source software: Implementation and management (Elsevier, 2004).\n14. Vigmond, E., Dos Santos, R. W., Prassl, A., Deo, M. & Plank, G. Solvers for the cardiac bidomain equations. Prog.\nbiophysics molecular biology 96, 3–18 (2008).\n15. Vázquez, M. et al. Alya: Multiphysics engineering simulation toward exascale. J. computational science 14, 15–27 (2016).\n16. Alnæs, M. et al. The FEniCS project version 1.5. Arch. numerical software 3 (2015).\n17. Plank, G. et al. The openCARP simulation environment for cardiac electrophysiology.Comput. Methods Programs Biomed.\n208, 106223 (2021).\n18. Africa, P. C. et al. lifex-ep: a robust and efficient software for cardiac electrophysiology simulations. BMC bioinformatics\n24, 389 (2023).\n19. Sachetto Oliveira, R. et al. Performance evaluation of GPU parallelization, space-time adaptive algorithms, and their\ncombination for simulating cardiac electrophysiology. Int. J. for Numer. Methods Biomed. Eng. 34, e2913 (2018).\n20. Gouvêa de Barros, B., Sachetto Oliveira, R., Meira Jr., W., Lobosco, M. & Weber dos Santos, R. Simulations of complex\nand microscopic models of cardiac electrophysiology powered by Multi-GPU platforms. Comput. Math. Methods Medicine\n2012, 824569 (2012).\n21. Sachetto Oliveira, R. et al. Comparing CUDA, OpenCL and OpenGL implementations of the cardiac monodomain\nequations. In Wyrzykowski, R., Dongarra, J., Karczewski, K. & Wa´sniewski, J. (eds.) Parallel Processing and Applied\nMathematics, 111–120 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).\n22. Amorim, R. M. & Weber dos Santos, R. Solving the cardiac bidomain equations using graphics processing units. J.\nComput. Sci. 4, 370–376 (2013).\n23. Kaboudian, A., Cherry, E. M. & Fenton, F. H. Large-scale interactive numerical experiments of chaos, solitons and fractals\nin real time via GPU in a web browser. Chaos, Solitons & Fractals 121, 6–29 (2019).\n24. Biasi, N., Seghetti, P., Parollo, M., Zucchelli, G. & Tognetti, A. A matlab toolbox for cardiac electrophysiology simulations\non patient-specific geometries. Comput. Biol. Medicine 185, 109529 (2025).\n25. Vergara, C. et al. A coupled 3D-1D numerical monodomain solver for cardiac electrical activation in the myocardium with\ndetailed Purkinje network. J. Comput. Phys. 308, 218–238 (2016).\n26. Sundnes, J. et al. On the computational complexity of the bidomain and the monodomain models of electrophysiology.\nAnnals biomedical engineering 34, 1088–1097 (2006).\n27. Mazumder, S. Numerical methods for partial differential equations: finite difference and finite volume methods (Academic\nPress, 2015).\n28. Gomes, J. M. et al. Uniformization method for solving cardiac electrophysiology models based on the markov-chain\nformulation. IEEE Transactions on Biomed. Eng. 62, 600–608 (2015).\n29. Zipes, D. P., Jalife, J. & Stevenson, W. G. Cardiac electrophysiology: from cell to bedside E-book (Elsevier Health\nSciences, 2017).\n30. Behradfar, E., Nygren, A. & Vigmond, E. J. The role of Purkinje-myocardial coupling during ventricular arrhythmia: a\nmodeling study. PLoS One 9, e88000 (2014).\n31. Wiedmann, R. T., Tan, R. C. & Joyner, R. W. Discontinuous conduction at Purkinje-ventricular muscle junction.Am. J.\nPhysiol. Circ. Physiol. 271, H1507–H1516 (1996).\n32. dos Santos, R. W. et al. ATX-II effects on the apparent location of M cells in a computational model of a human left\nventricular wedge. J. cardiovascular electrophysiology 17, S86–S95 (2006).\n33. Vigmond, E. J. & Clements, C. Construction of a computer model to investigate sawtooth effects in the Purkinje system.\nIEEE transactions on biomedical engineering 54, 389–399 (2007).\n34. Bishop, M. J. & Plank, G. Bidomain ECG simulations using an augmented monodomain model for the cardiac source.\nIEEE transactions on biomedical engineering 58, 2297–2307 (2011).\n35. Niederer, S. A. et al. Verification of cardiac tissue electrophysiology simulators using an N-version benchmark.Philos.\nTransactions Royal Soc. A 369, 4331–4351 (2011).\n11/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\n36. ten Tusscher, K. H., Noble, D., Noble, P.-J. & Panfilov, A. V . A model for human ventricular tissue.Am. J. Physiol. Circ.\nPhysiol. 286, H1573–H1589 (2004).\n37. Tomek, J. et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease,\nand drug block. Elife 8, e48890 (2019).\n38. Rush, S. & Larsen, H. A practical algorithm for solving dynamic membrane equations. IEEE Transactions on Biomed.\nEng. 389–392 (1978).\n39. Banerjee, A. et al. A completely automated pipeline for 3D reconstruction of human heart from 2D cine magnetic resonance\nslices. Philos. Transactions Royal Soc. A 379, 20200257 (2021).\n40. Camps, J. et al. Inference of ventricular activation properties from non-invasive electrocardiography. Med. Image Analysis\n73, 102143 (2021).\n41. Camps, J. et al. Digital Twinning of the Human Ventricular Activation Sequence to Clinical 12-lead ECGs and Magnetic\nResonance Imaging Using Realistic Purkinje Networks for in Silico Clinical Trials. Med. Image Analysis 94, 103108\n(2024).\n42. Camps, J. et al. Harnessing 12-lead ECG and MRI data to personalise repolarisation profiles in cardiac digital twin models\nfor enhanced virtual drug testing. Med. Image Analysis 100, 103361 (2025).\n43. Doste, R., Coppini, R. & Bueno-Orovio, A. Remodelling of potassium currents underlies arrhythmic action potential\nprolongation under beta-adrenergic stimulation in hypertrophic cardiomyopathy. J. Mol. Cell. Cardiol. 172, 120–131\n(2022).\n44. Trovato, C. et al. Human Purkinje in silico model enables mechanistic investigations into automaticity and pro-arrhythmic\nabnormalities. J. Mol. Cell. Cardiol. 142, 24–38 (2020).\n45. Li, J., Zhang, H. & Boyett, M. Numerical analysis of conduction of the action potential across the Purkinje fibre-ventricular\nmuscle junction. In Computing in Cardiology, 265–268 (2016).\n46. Mendez, C., Mueller, W. J. & Urguiaga, X. Propagation of impulses across the Purkinje fiber-muscle junctions in the dog\nheart. Circ. research 26, 135–150 (1970).\n47. Matsuda, K. Configuration of the transmembrane potential of the Purkinje-ventricular fiber junction and its analysis.\nElectrophysiol. Ultrastruct. Hear. 177–187 (1967).\nData Availability\nThe open source cardiac solver is publicly available at https://github.com/rsachetto/MonoAlg3D_C. All the\nnecessary configuration files, custom functions, post-processing scripts used during the current study, as well as, the biventricular\nmesh and Purkinje networks used for the cardiac digital twin application will be publicly available at a Zenodo repository to\nallow reproducibility once the manuscript has been accepted.\nAcknowledgements\nThis work was funded by a Wellcome Trust fellowship in Basic Biomedical Sciences to B.R. (214290/Z/18/Z), the EPSRC\nproject CompBioMedX (EP/X019446/1), the CompBioMed2 Centre of Excellence in Computational Biomedicine grant\nagreements No. 675451 and No. 823712, R.S.O. acknowledges support from Fapemig grant No. APQ-00748-18 and UFSJ.\nR.W.S. acknowledges support from Fapemig grant APQ-02445-24 and UFJF. A.B.O. acknowledges support from UK Research\nand Innovation grant No. 10110728. The U.S. Department of Energy’s (DOE) Innovative and Novel Computational Impact\non Theory and Experiment (INCITE) Program awarded access to Polaris, under contract No. DE-AC02-06CH11357. For the\npurpose of open access, the authors have applied a Creative Commons Attribution (CC BY) public copyright licence to any\nAuthor Accepted Manuscript version arising from this submission.\nAuthor contributions statement\nL.A.B.: Conceptualization, Methodology, Software, Investigation, Formal analysis, Writing – original draft, Writing – review\n& editing. R.S.O.: Conceptualization, Methodology, Software, Investigation, Formal analysis, Writing – original draft, Writing\n– review & editing. J.C.: Software, Formal analysis, Writing – review & editing. Z.J.W.: Software, Formal analysis, Writing –\nreview & editing. R.D.: Data curation, Formal analysis, Writing – review & editing. A.B.O.: Formal analysis, Writing – review\n& editing. R.W.S.: Conceptualization, Methodology, Investigation, Formal analysis, Writing – original draft, Writing – review\n& editing. B.R.: Conceptualization, Formal analysis, Funding acquisition, Resources, Writing – review & editing.\n12/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint \n\nAdditional information\nSupplementary material accompanies this paper.\nCompeting interests: The authors declare no competing interests.\n13/13\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 9, 2025. ; https://doi.org/10.1101/2025.04.09.647733doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}