Digital twinning of cardiac electrophysiology for congenital heart disease

2.1. Pre-processing

Figure 1 (first row) shows the heart–torso model of a 7-year-old female paediatric patient with HLHS constructed from the computerized tomography (CT) scan of the patient using our semi-automatic model construction pipeline [29]. This patient has a severely underdeveloped (hypoplastic) left ventricle (LV) and a dominant right ventricle (RV), which is connected to the aorta. Figure 2 shows the CHD captured in this HLHS patient, compared with the anatomy of a healthy patient.

• Download figure
• Open in new tab
• Download PowerPoint

2.2. Cardiac electrophysiology modelling

We run 200 numerical simulations on the patient-specific heart–torso geometry (see figure 1, second row), spanning seven relevant electrophysiology parameters of the physics-based model at the microscopic scale and organ level. We collect the corresponding in silico 12-lead ECGs.

In table 1, we report descriptions, ranges, units and references for the seven model parameters that we explore via Latin hypercube sampling for the dataset generation. For ionic conductances and myocardial conductivities, we consider a healthy baseline value from the literature and we define an interval by multiplying this value by 0.5 (lower bound) and 2.0 (upper bound). In this way, we cover a wide range that includes both healthy and pathological conditions. Similarly, for the Purkinje conductivities, we use an interval spanning both healthy cases and different cardiovascular diseases [38]. For the stimulation time of the left Purkinje network, we employ a broad range for possible ventricular dyssynchrony in patients with single ventricle physiology [27].

In figure 3, we depict the ensemble of the resulting in silico 12-lead ECGs together with the clinical recordings. We point out that the patient-specific 12-lead ECGs are contained within the pseudopotential variability spanned by the electrophysiology simulations, manifesting various morphologies in the QRS complex, that is ventricular depolarization, and T wave, that is ventricular repolarization. The patient diagnosis reports rhythm disorders, atrial enlargement, left and right ventricular hypertrophy, along with severe abnormalities in the ECGs. Specifically, there are signs of prolonged PR interval, ST segment depression and T wave inversion.

• Download figure
• Open in new tab
• Download PowerPoint

In figure 4, we show the simulated spatio-temporal transmembrane potential evolution on the patient-specific paediatric model for a single instance of model parameters. Specifically, we always simulate the sinus rhythm behavior over a cardiac cycle. Figure 4 focuses on the ventricular depolarization phase, where the electric signal propagates from the one-dimensional Purkinje network at the two endocardia towards the myocardium, as well as the ventricular repolarization phase, when the transmembrane potential comes back to its resting state (i.e. approximately −90 mV).

• Download figure
• Open in new tab
• Download PowerPoint

2.3. Branched latent neural maps

We train BLNMs, which are represented by feed-forward partially connected NNs, to encode the temporal dynamics of the 12-lead pseudo-ECGs computed with the physics-based mathematical model while also covering model variability from the cellular to the tissue level (see figure 1, third row). Once trained, BLNMs act as a surrogate model for cardiac electrophysiology function that can be queried on new parameter instances (within the training range) to provide faster than real-time in silico 12-lead ECGs.

In order to identify the optimal set of BLNM hyperparameters, which are the number of layers, number of neurons, number of states and disentanglement level in the NN structure, we employ a ${\displaystyle {\displaystyle K}}$-fold ($K=5$) cross-validation over 150 multi-scale physics-based electrophysiology simulations. The hyperparameter search space is given by a four-dimensional hypercube, where we run 50 instances of Latin hypercube sampling and we pick the BLNM configuration providing the lowest generalization error. For each configuration of hyperparameters, we sample the dataset with a fixed time step of $\mathrm{\Delta }t=5$ ms, and we perform 10 000 iterations of the second-order Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimizer. In table 2, we report the initial hyperparameter ranges for tuning and the final optimized values.

Then, we train a final optimized BLNM encompassing the whole dataset of 150 multi-scale physics-based electrophysiology simulations using 50 000 BFGS iterations. The non-dimensional MSE on a testing set comprised of 50 additional numerical simulations unseen during the training stage, and Latin hypercube sampled from the same parameter space in table 1, is equal to ${\displaystyle {\displaystyle 5\times {10}^{-4}}}$.

2.4. Parameter estimation

We employ the optimized BLNM to find an initial guess for the seven model parameters that results in a computational pseudo-ECG that best matches the clinically observed 12-lead ECG dynamics of the CHD patient. The estimated model parameters are reported in table 3.

2.5. Sensitivity analysis

Starting from the parameter calibration shown in table 3, we compute Shapley effects for each model parameter for cardiac electrophysiology, assuming independence among them as they act in different terms and equations of the physics-based mathematical model (see figure 1, third row). In figure 5, we show how each parameter contributes to matching electrophysiology simulations with the clinical 12-lead ECGs, that is, in the minimization of the MSE between BLNM outputs and our patient-specific observations. The sodium current conductance $G_{\mathrm{Na}}$ plays a dominant role, followed by the L-type calcium ion channel conductance $G_{\mathrm{CaL}}$ and the different conductivities $D_{\mathrm{ani}}$, $D_{\mathrm{iso}}$ and $D_{\mathrm{purk}}$. Noteworthy, the interventricular activation dyssynchrony $t_{\mathrm{LV}}^{\mathrm{stim}}$ plays a minor role in the calibration process. This is motivated by the dimension of the right ventricle, which mostly dictates the activation sequence with respect to the small (underdeveloped) left ventricle.

• Download figure
• Open in new tab
• Download PowerPoint

2.6. Uncertainty quantification

In figures 6 and 7, we show the results of our inverse uncertainty quantification, where we quantify how uncertainty in matching 12-lead ECGs propagates to uncertainty in the estimated model parameters. We account for both BLNM surrogate modelling error, via GPs, and measurement error in the clinical recordings.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

From figure 6, we see that the posterior distributions of all model parameters ${\displaystyle {\displaystyle {\mathit{\theta }}_{\mathrm{E}\mathrm{P}}}}$, along with the correlation length $l_{\mathrm{GP}}$ and amplitude ${\sigma }_{\mathrm{GP}}$, converged well towards highly similar unimodal distributions across all chains. The average value of ${\sigma }_{\mathrm{GP}}^2$ is approximately equal to 0.08, which is two orders of magnitude higher than the BLNM testing error (${\displaystyle {\displaystyle 5\times {10}^{-4}}}$), as the GP encodes the maximum BLNM uncertainty from each single lead individually and by also considering all possible correlations among the 12 leads, given the full covariance matrix in the multi-variate normal distribution (see §3.6).

In figure 7, we depict the clinical versus in silico 12-lead ECGs, generated with the BLNM over the posterior distribution of model parameters. We see that the numerical simulations are in good agreement with the patient-specific recordings and show small variability between the 5 s.d. from the average value.

2.7. In silico scenarios

We consider different in silico scenarios of clinical interest using our personalized computational model. Specifically, we run 3 one-dimensional–three-dimensional electrophysiology simulations under three different conditions: first, we use the calibrated model parameters ${\displaystyle {\displaystyle {\mathit{\theta }}_{\mathrm{E}\mathrm{P}}}}$, that is, the peak of the posterior distribution in figure 6, then we employ the same computational model under the assumption of either left bundle branch block or right bundle branch block, that is, by removing the left (respectively, right) Purkinje fascicles. In figure 8, we show the results of these three numerical simulations by depicting the activation and repolarization maps for this HLHS paediatric case. We notice that, for this patient-specific geometry, the role of the Purkinje network in the LV is very limited and the activation sequence is highly similar with and without a full left bundle branch block. This means that the left Purkinje fascicle in this patient is either absent or present but fully inhibited or severely delayed, with respect to the right Purkinje fascicle.

• Download figure
• Open in new tab
• Download PowerPoint

2.8. Computational costs

In table 4, we detail the computational costs and resources required by each step of the digital twinning process. The most expensive part resides in the physics-based computational electrophysiology modelling dataset generation, which makes use of high-performance computing given the stiffness and complexity of the underlying mathematical model. On the other hand, training a BLNM and employing it for robust Bayesian parameter estimation and sensitivity analysis are more tractable tasks that can be carried out within a few hours or minutes on a local machine. Using the physics-based mathematical model throughout parameter calibration with uncertainty quantification and sensitivity analysis would be computationally intractable and unaffordable given the extensive number of queries that Shapley effects and HMC require to show robustness and convergence in the provided results (see §3).

3.1. Pre-processing

We reconstruct the heart–torso model from the CT scan of a 7-year-old female paediatric patient with HLHS in a semi-automatic manner [29]. Namely, we train a NN based on the classic UNet architecture [41] to automatically segment the myocardium from CT images of patients affected by congenital heart disease. The UNet is trained using a publicly available dataset [42] that provided CT images and ground-truth segmentation for 110 patients with age between 1 month and 40 years, combined with our private HLHS dataset containing the images and segmentation of five patients. Given the intrinsic segmentation challenges of cardiac structures in both young and CHD patients [43], we subsequently examine and improve the UNet-produced segmentations to more closely match with the CT scan. We automatically extract the surface meshes from the segmentations using the marching cube algorithm [44] and truncate the base myocardium above a manually identified plane to create a biventricular surface mesh. We choose this plane so that it is normal to the long axis of the ventricles and located below the aortic and atrioventricular annulus. We manually adjust the location of the plane to completely remove the irregular annulus geometries while keeping most of the ventricles (see figure 2). We note that this choice is common to other electrophysiological studies [22,45] and facilitates the generation of the fibre field by rule-based algorithms [46]. We subsequently use TetGen [47] to create the tetrahedral volumetric mesh with a maximum edge size of 1 mm [27,32]. The torso model is created semi-automatically from the images using threshold- and region-growing-based segmentation methods.

3.2. Cardiac electrophysiology modelling

We detail the physics-based mathematical model, along with its numerical discretization, that is employed to perform electrophysiology simulations on the HLHS paediatric patient.

3.3. Branched latent neural maps

We construct a geometry-specific surrogate model of cardiac function by building a feed-forward partially connected NN that explores the variability of our physics-based electrophysiology model detailed in §3.2 while structurally separating the role of temporal $t$ and functional ${\displaystyle {\displaystyle {\mathit{\theta }}_{\mathrm{E}\mathrm{P}}}}$ parameters. This recent scientific machine-learning tool, proposed by Salvador & Marsden [32], allows for different levels of disentanglement between inputs and outputs. The surrogate model reads

Weights and biases ${\displaystyle {\displaystyle \mathbf{w}\in {\mathbb{R}}^{N_{\mathrm{w}}}}}$ encode the algebraic structures of a feed-forward partially connected NN, which represents a map ${\displaystyle {\displaystyle \mathcal{B}\mathcal{L}\mathcal{N}\mathcal{M}:{\mathbb{R}}^{1+N_{\mathcal{P}}}\to {\mathbb{R}}^{N_{\mathrm{z}}}}}$ from time $t$ and ${\displaystyle {\displaystyle N_{\mathcal{P}}=7}}$ cell-to-organ scale electrophysiology parameters ${\displaystyle {\displaystyle {\mathit{\theta }}_{\mathrm{E}\mathrm{P}}=[G_{\mathrm{C}\mathrm{a}\mathrm{L}},G_{\mathrm{N}\mathrm{a}},G_{\mathrm{K}\mathrm{r}},D_{\mathrm{a}\mathrm{n}\mathrm{i}},D_{\mathrm{i}\mathrm{s}\mathrm{o}},D_{\mathrm{p}\mathrm{u}\mathrm{r}\mathrm{k}},t_{\mathrm{L}\mathrm{V}}^{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}{]}^T\in \mathbf{\Theta }\subset {\mathbb{R}}^{N_{\mathcal{P}}}}}$ to an output vector ${\displaystyle {\displaystyle \mathbf{z}(t)=[{\mathbf{z}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{s}}(t),{\mathbf{z}}_{\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}}(t){]}^T\in {\mathbb{R}}^{N_{\mathrm{z}}}}}$. This vector contains in silico precordial and limb leads recordings ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{s}}(t)=[V_1(t),V_2(t),V_3(t),V_4(t),V_5(t),V_6(t),LA(t),RA(t),F(t){]}^T\in {\mathbb{R}}^9}}$, where we use the original $LA(t)$, $RA(t)$ and $F(t)$ limb leads in place of the bipolar and augmented limb leads in order to reduce the dimensionality of the output. Indeed, we reconstruct $I(t)$, $II(t)$, $III(t)$, $aVL(t)$, $aVR(t)$ and $aVF(t)$ a posteriori by means of equations (3.3) and (3.4). Furthermore, vector ${\displaystyle {\displaystyle \mathbf{z}(t)}}$ leverages some ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}}(t)}}$ latent variables that enhance the learned temporal dynamics by acting in regions with steep gradients [32].

We perform nonlinear optimization with the BFGS algorithm to tune NN parameters. In particular, we monitor the MSE of surrogate versus physics-based ECG pseudopotentials to find an optimal set of weights and biases ${\displaystyle {\displaystyle \mathbf{w}}}$, that is,

where ${\displaystyle {\displaystyle {\tilde{\mathbf{z}}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{s}}(t)\in [-1,1{]}^9}}$ represents BLNM outputs and ${\displaystyle {\displaystyle {\tilde{\mathbf{z}}}_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}(t)\in [-1,1{]}^9}}$ defines the physics-based numerical simulations, both in non-dimensional form. Time $\tilde{t}\in [0,1]$ and model parameters ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}\in [-1,1{]}^{N_{\mathcal{P}}}}}$ are also normalized during the training and testing phases. We refer to Salvador & Marsden [32] for a detailed description of all the properties related to BLNMs that enable them to effectively learn complex physical processes.

3.4. Parameter estimation

We employ our trained BLNM to find a set of model parameters ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}}}$ that matches ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{E}\mathrm{C}\mathrm{G}}(t)\in {\mathbb{R}}^{12}}}$ with ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}(t)\in {\mathbb{R}}^{12}}}$. Here, ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{E}\mathrm{C}\mathrm{G}}(t)}}$ is the vector of BLNM physical outputs ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{s}}(t)}}$ manipulated according to equations (3.3) and (3.4) to generate the full 12-lead ECGs, and ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}(t)}}$ is the clinically measured ECG data vector. In particular, we perform derivative-free optimization by employing the Nelder–Mead method [54], where we specify a loss function given by the MSE of the mismatch between the trained surrogate versus clinical ECG potentials, that is,

which leads to a set of calibrated model parameters ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}^{\mathrm{N}\mathrm{M}}}}$ and corresponding 12-lead ECGs ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{E}\mathrm{C}\mathrm{G}}^{\mathrm{N}\mathrm{M}}(t)}}$. We initialize our optimization algorithm with a random set of model parameters ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\in [-1,1{]}^{N_{\mathcal{P}}}}}$. We repeat the optimization process 100 times and we average the model parameters obtained during the different trials in order to get ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}^{\mathrm{N}\mathrm{M}}}}$.

3.5. Sensitivity analysis

We perform a variance-based sensitivity analysis using Shapley effects [33] in order to quantify the importance of each model parameter in fitting patient-specific 12-lead ECGs during the inference process.

Specifically, we employ Sklar’s theorem [55] to define the input multi-variate distribution, which is given by a Gaussian copula and a series of ${\displaystyle {\displaystyle N_{\mathcal{P}}}}$ marginals defined by standard normal distributions centred in ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}^{\mathrm{N}\mathrm{M}}}}$, that is, ${\displaystyle {\displaystyle \mathcal{N}\left({\tilde{\theta }}_{\mathrm{E}\mathrm{P}}^{\mathrm{N}\mathrm{M},\mathrm{i}},0.2\right)\text{for}i=1,...,N_{\mathcal{P}}}}$.

Due to the high computational costs associated with testing all the different combinations of the features, we consider the random (rather than the exact) version of the algorithm to compute Shapley values. We monitor the expected marginal contribution of each model parameter to the BLNM prediction with respect to observations, that is, the MSE of equation (3.7). We fix 2000 permutations, 500 bootstrapped samples and 50 samples to estimate conditional variance for three times.

3.6. Uncertainty quantification

We employ a BLNM within HMC [34] to calibrate model parameters and to perform inverse uncertainty by matching observed 12-lead ECGs from patient-specific recordings. HMC is a Markov chain Monte Carlo (MCMC) method that aims at finding an approximation of the posterior distribution ${\displaystyle {\displaystyle \mathbb{P}({\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}|\mathbf{x})}}$, given a certain prior probability distribution ${\displaystyle {\displaystyle \mathbb{P}({\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}})}}$ with respect to the model parameters in non-dimensional form ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}\in [-1,1{]}^{N_{\mathcal{P}}}}}$. Specifically, we employ the No-U-Turn Sampler (NUTS) extension of HMC, which automatically adapts the number of steps to estimate the posterior distribution [35]. This algorithm, which shares and enhances some of the features of sequential [56] and differential evolution [57] MCMC, works well with high-dimensional target distributions, possibly presenting correlated dimensions. Moreover, HMC reaches convergence using a reduced amount of samples with respect to vanilla MCMC [35]. For further details about the mathematical derivation of HMC and its application to cardiac simulations, we refer to Salvador et al. [20].

We run four chains with 1000 adaptation samples in the warm-up phase and 1000 effective samples to estimate the posterior distribution, with a fixed $90\%$ acceptance rate. For all model parameters, we consider prior distributions

where ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}^{\mathrm{N}\mathrm{M}}\in [-1,1{]}^{N_{\mathcal{P}}}}}$ is the initial guess obtained with the Nelder–Mead method. We always make sure that model parameters reside within the $[-1,1]$ range. We set $\iota =0.2$. Even though NUTS allows for many different initialization protocols, such as maximum a posteriori (MAP) or maximum likelihood estimation (MLE), we consider an initial random seed for each chain. This is motivated by the sensitivity of MAP and MLE over multiple runs, especially when several model parameters are calibrated with respect to noisy or highly varying time-dependent QoIs, which is the case for ECG recordings.

Several sources of uncertainty can be considered. These include model uncertainties (e.g. the discrepancy between the actual physical phenomenon and the high-fidelity model), the discretization error introduced when solving the differential equations, the surrogate modelling error of reduced-order models and the measurement errors that intrinsically affect clinical ECG recordings (i.e. the sensitivity of the instrument used during the clinical test, variations in lead placement position by clinicians and patient-specific factors such as breathing and motion). However, in our inverse uncertainty quantification process, we only include the measurement error and the approximation error introduced by the transition from the high-fidelity model to the BLNM-based surrogate model. In particular, we consider a multi-variate normal distribution centred in the BLNM predictions ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{E}\mathrm{C}\mathrm{G}}(t)}}$ for the given patient-specific observations ${\displaystyle {\displaystyle {\mathbf{z}}_{\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}(t)}}$, which reads

${\sigma }_{\mathrm{meas}}=0.1$ is the a priori fixed standard deviation dictating the measurement error [58,59], whereas ${\displaystyle {\displaystyle k(\stackrel{\mathbf{~}}{\mathit{t}},{\stackrel{\mathbf{~}}{\mathit{t}}}^{\mathrm{\prime }};{\sigma }_{\mathrm{G}\mathrm{P}},l_{\mathrm{G}\mathrm{P}})={\sigma }_{\mathrm{G}\mathrm{P}}^2\exp \left(\frac{-||\stackrel{\mathbf{~}}{\mathit{t}}-{\stackrel{\mathbf{~}}{\mathit{t}}}^{\mathrm{\prime }}|{|}^2}{2l_{\mathrm{G}\mathrm{P}}^2}\right)}}$ is the exponentiated quadratic kernel of a zero-mean GP ${\displaystyle {\displaystyle \mathcal{G}\mathcal{P}(\mathbf{0},k(\stackrel{\mathbf{~}}{\mathit{t}},{\stackrel{\mathbf{~}}{\mathit{t}}}^{\mathrm{\prime }};{\sigma }_{\mathrm{G}\mathrm{P}},l_{\mathrm{G}\mathrm{P}}))}}$ [60]. Amplitude ${\sigma }_{\mathrm{GP}}\sim (0.01,1.0)$ and correlation length $l_{\mathrm{GP}}\sim (0.01,1.0)$ are additional hyperpriors tuned during HMC to quantify the surrogate modelling error, which may change according to the specific observation. Vectors ${\displaystyle {\displaystyle \stackrel{\mathbf{~}}{\mathit{t}}}}$ and ${\displaystyle {\displaystyle {\stackrel{\mathbf{~}}{\mathit{t}}}^{\mathrm{\prime }}}}$ represent discrete time points in the $[0,1]$ interval.

A full covariance matrix in the multi-variate normal distribution allows us to model the correlation among different leads. We evaluate the convergence of the HMC chains by checking that the Gelman–Rubin diagnostic provides a value less than 1.1 for all the model parameters ${\displaystyle {\displaystyle {\tilde{\mathit{\theta }}}_{\mathrm{E}\mathrm{P}}}}$, $l_{\mathrm{GP}}$ and ${\sigma }_{\mathrm{GP}}$ [61–63].

3.7. Software and hardware

All electrophysiology simulations are performed at the Stanford Research Computing Center using svFSIplus [64], a C++ high-performance computing multi-physics and multi-scale finite element solver for cardiac and cardiovascular modelling. This solver is part of the SimVascular software suite for patient-specific cardiovascular modelling [65].

We train the NNs by using BLNM.jl [32], Julia library for scientific machine learning which is publicly available under MIT License at https://github.com/StanfordCBCL/BLNM.jl. This library leverages Hyperopt.jl [66] for parallel hyperparameter optimization by combining the Message Passing Interface with Open Multi-Processing (OpenMP) on physical and virtual cores, respectively.

We perform sensitivity analysis and parameter estimation with uncertainty quantification using GlobalSensitivity.jl [67] and Turing.jl [68], respectively, which both exploit OpenMP and vectorized operations to speed-up computations. The code for sensitivity analysis and Bayesian parameter estimation is available within BLNM.jl as a test case.

Furthermore, this public repository contains the dataset encompassing all the electrophysiology simulations used for the training and testing phases, along with the patient-specific 12-lead ECGs.