Digital twinning of cardiac electrophysiology for congenital heart disease

In recent years, blending mechanistic knowledge with machine learning has had a major impact in digital healthcare. In this work, we introduce a computational pipeline to build certified digital replicas of cardiac electrophysiology in pediatric patients with congenital heart disease. We construct the patient-specific geometry by means of semi-automatic segmentation and meshing tools. We generate a dataset of electrophysiology simulations covering cell-to-organ level model parameters and utilizing rigorous mathematical models based on differential equations. We previously proposed Branched Latent Neural Maps (BLNMs) as an accurate and efficient means to recapitulate complex physical processes in a neural network. Here, we employ BLNMs to encode the parametrized temporal dynamics of in silico 12-lead electrocardiograms (ECGs). BLNMs act as a geometry-specific surrogate model of cardiac function for fast and robust parameter estimation to match clinical ECGs in pediatric patients. Identifiability and trustworthiness of calibrated model parameters are assessed by sensitivity analysis and uncertainty quantification.


Introduction
The combination of physics-based and statistical modeling in cardiovascular medicine has the potential to shape the future of cardiology [9].In this framework, a synergistic use of multiphysics and multiscale mathematical models for cardiac function [13,18,39,41,57] and machine learningbased methods, such as Gaussian processes emulators [30,37,59] and Neural Networks (NNs) [35,47], enables the design of efficient computational tools that are compatible with the computer resources and time frames required in clinical applications.In the foreseeable future, a continuous, bi-directional interaction between patient-specific data and Artificial Intelligence-enriched computer models incorporating biophysically detailed and anatomically accurate knowledge would enable the vision of precision medicine [34,38,44].Personalized treatment and surgical planning may be delivered by leveraging different mathematical methods, such as sensitivity analysis, parameter inference and uncertainty quantification [19,25,49,54].
Several mathematical tools have been proposed to better understand and treat different groups of adult cardiac pathologies [34].Electrophysiology simulations play an important role in the assessment of rhythm disorders.They are used for cardiac resynchronization therapy [58], arrhythmia risk stratification [2,50], and definition of optimal ablation strategies [43].Nevertheless, in silico numerical simulations and treatment modalities in pediatrics and congenital heart disease are less common or even not established [8,24,62].
Congenital heart defects (CHDs) are the most common birth defects and are characterized by cardiac anatomical abnormalities that can severely impact cardiac function [29].Patients with CHDs often have a unique and peculiar combination of cardiac defects that warrant personalized treatment planning in clinically-relevant time frames.Digital twinning of cardiac function thus holds particular promise for these patients [62].
In this work we introduce a novel digital twin of a pediatric patient with hypoplastic left heart syndrome (HLHS), a complex form of CHD where the left ventricle of the patient is severely underdeveloped, leading to a number of morbidities and elevated mortality risk.Our pipeline seamlessly combines: • Semi-automatic segmentation and mesh generation tools suited for pediatric patients with CHD [27] • A physiologically-based mathematical formulation of cardiac electrophysiology deriving from the monodomain equation [44,51] coupled with the ten Tusscher-Panfilov [63] ionic model • A recently proposed scientific machine learning method, namely Branched Latent Neural Maps (BLNMs) [53], to build an accurate and efficient dynamic surrogate model of cardiac function • Shapley effects [56] and Hamiltonian Monte Carlo (HMC) [5,21] to perform patient-specific sensitivity analysis, fast and robust parameter estimation with uncertainty quantification while matching clinical 12-lead electrocardiograms (ECGs) This digital twin can be employed to simulate different scenarios of clinical interest in silico, as HLHS patients may experience different forms of electrophysiological comorbidities [14], including ventricular arrhythmias and dyssynchrony [62].Therefore, personalized electrophysiology simulations may provide virtual pre-and post-operative guidance in this understudied patient cohort [32].

Results
In Figure 1 we depict our computational pipeline to build digital twins of cardiac electrophysiology for congenital heart disease patients.This pipeline covers all the relevant aspects of digital twinning: image segmentation and mesh generation, mathematical and numerical physics-based modeling, surrogate model training, sensitivity analysis and robust parameter calibration with uncertainty quantification.

Pre-processing
Figure 1 (first row) shows the heart-torso model of a 7-year-old female pediatric patient with HLHS constructed from the computerized tomography (CT) scan of the patient using our semi-automatic model construction pipeline [27].

Segmentation Surface Processing Meshing
Higher input-output disentanglement Lower input-output disentanglement  We reconstruct the patient-specific geometry with HLHS from imaging.We generate a dataset of electrophysiology simulations encompassing cell-toorgan variability in model parameters.We train a BLNM that effectively reproduces 12-lead ECGs while covering model variability.We employ the BLNM for digital twinning.

Cardiac electrophysiology modeling
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, and units for the seven model parameters that we explore via latin hypercube sampling for the dataset generation.
In Figure 2 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 pseudopotentials 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.
In Figure 3 we show the simulated spatio-temporal transmembrane potential evolution on the patient-specific pediatric model for a single instance of model parameters.Specifically, we always simulate the sinus rhythm behavior over a cardiac cycle.Figure 3 focuses on the ventricular depolarization phase, where the electric signal propagates from the 1D 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).We train BLNMs, which are represented by feedforward 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 to provide faster than real-time in silico 12-lead ECGs.

Branched Latent Neural Maps
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 Kfold (K = 5) cross validation over 150 multiscale 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 ∆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 hyperparameters ranges for tuning and the final optimized values.Then, we train a final optimized BLNM encompassing the whole dataset of 150 multiscale physics-based electrophysiology simulations using 50,000 BFGS iterations.The non-dimensional Mean Square Error (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 5 • 10 −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.
Even though the relative heart orientation and lead placements significantly influence ECGs [66], and as such may require additional parameters to calibrate, the information retrieved from the CT scan and patient diagnosis allow us to determine these quantities with a small degree of uncertainty.This motivates our focus in the estimation process on the cell-to-organ level electrophysiology model parameters which are assessed as important in previous studies [10,60].

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 4 we show how each parameter contributes in matching electrophysiology simulations with the clinical 12-lead ECGs, i.e. in the minimization of the MSE between BLNM outputs and our patient-specific observations.The sodium current conductance G Na plays a dominant role, followed by the L-type calcium ion channel conductance G CaL and the different conductivities D ani , D iso and D purk .Noteworthy, the interventricular activation dyssynchrony t stim LV 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.

Uncertainty quantification
In Figures 5 and 6 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 modeling error, via Gaussian processes (GPs), and measurement error in the clinical recordings.
From Figure 5, we see that the posterior distributions of all model parameters θ EP , along with the correlation length l GP and amplitude σ GP , converged well towards highly similar unimodal distributions across all chains.The average value of σ 2 GP is approximately equal to 0.08, which is two orders of magnitude higher than the BLNM testing error (5 • 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 multivariate normal distribution (see Section 3.6).
In Figure 6 we depict the clinical vs. 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 five standard deviations from the average value.

In silico clinical trial
In Figure 7 we show the results of three numerical simulations in which we analyze different scenarios of clinical interest pre-operatively on the HLHS pediatric patient.Specifically, we depict activation and repolarization times for the electrophysiology simulation with the calibrated model parameters θ EP and by inducing either a left or a right bundle branch block, where the left (respectively, right) Purkinje fascicles are inhibited.We notice that, for this patient-specific case, the role of the Purkinje network in the left ventricle is very limited and that the activation sequence is highly similar with and without a full left bundle branch block.

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 modeling 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 Figure 6: Inverse uncertainty quantification: matching clinical data.Clinical recordings (dashed, black) and mean estimation (red, solid) of 12-lead ECGs for the HLHS pediatric patient via HMC.Light red encompasses the variability between mean minus/plus five standard deviations.tasks that can be carried out within a few hours or minutes on a local machine.Using the physicsbased 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 Methods section).

Pre-processing
We reconstruct the heart-torso model from the CT scan of a 7-year-old female pediatric patient with HLHS in a semi-automatic manner [27].Namely, we train a NN based on the classic UNet architecture [23] to automatically segment the myocardium from CT images.The UNet is trained using a publicly available dataset [67] 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 Table 4: Computational resources.Summary of the computational times and resources to generate the electrophysiology simulations with the physics-based model, to train the BLNM, to compute Shapley values for sensitivity analysis and to perform Bayesian parameter estimation with uncertainty quantification on 12-lead ECGs.
containing the images and segmentation of 5 patients.Given the intrinsic segmentation challenges of cardiac structures in both young and CHD patients [40], 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 [31] and truncate the base myocardium above a manually identified plane to create a biventricular surface mesh.We subsequently use TetGen [55] to create the tetrahedral volumetric mesh with a maximum edge size of 1 mm [53,62].The torso model is created semi-automatically from the images using thresholdand region-growing-based segmentation methods.Images and associated clinical data were obtained under an IRB-approved protocol at Stanford University.

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

Mathematical model
We consider the monodomain equation [44] coupled with the ten Tusscher-Panfilov ionic model [63] to describe the electric behavior in the heart-Purkinje system.This system of differential equations reads: We always simulate a single heartbeat by fixing a final time T = T HB = 600 ms.The computational domain Ω = Ω purk ∪ Ω myo is represented by the one-way coupled 1D Purkinje network and 3D biventricular patient-specific geometry.
Transmembrane potential Φ defines the electric signal at the Purkinje and myocardial level.The ten Tusscher-Panfilov ionic model is endowed with 18 variables, which are split in two different subsets.First, there is a vector w = (w 1 , . . ., w M ) (M = 12) of ion channel gating variables, which are probability density functions representing the fraction of open channels across the membrane of a single cardiac cell.Then, there is a vector z = (z 1 , . . ., z P ) (P = 6) of concentration variables representing relevant ionic species, such as sodium N a + , intracellular calcium Ca 2+ and potassium K + , which all play a major role in the metabolic processes [4], dictating heart rhythmicity or sarcomere contractility, and are generally targeted by pharmaceutical therapies.The specific mathematical formulation of the ten Tusscher-Panfilov ionic model defines the ordinary differential equations for H(Φ, w, z) and G(Φ, w, z), which describe the dynamics of gating variables and ionic concentrations respectively, along with the ionic current I ion (Φ, w, z) [63].An external applied current I app (x, t) fires the electric signal in the Purkinje fibers.
The diffusion tensor is expressed as , where f 0 introduces the biventricular fiber field [42,51].D ani , D iso , D purk ∈ R + dictate the anisotropic, isotropic and Purkinje conductivities, respectively.
The homogeneous Neumann boundary conditions prescribed at ∂Ω define the condition of an electrically isolated domain, where n is the outward unit normal vector to the boundary.
Following [51], the extracellular potential Φ e defining the ECGs is computed in each lead location x e as: where e = {V 1 , V 2 , V 3 , V 4 , V 5 , V 6 } and e = {LA, RA, F } define six precordial leads and three limb leads located on the pediatric patient-specific torso model, shown in colored and black dots in Figure 1 (second row), respectively.From these lead locations, we computationally reconstruct three bipolar limb leads as: and three augmented limb leads as: The resulting set ECG = {V 1 , V 2 , V 3 , V 4 , V 5 , V 6 , I, II, III, aV L, aV R, aV F } of computational pseudopotentials defines a comprehensive 12-lead ECG representation of the electrical activity in the patient-specific heart.

Numerical discretization
We employ linear Finite Elements to discretize the spatial domain Ω in Equation (3.2.2).The tetrahedral tessellation defining the biventricular mesh has 933,916 cells and 158,277 DOFs, with a maximum mesh size of h = 1 mm.The left and right Purkinje bundles within the ventricular endocardia are generated by employing the fractal tree and projection algorithm proposed in [48], starting from the atrioventricular node.These left and right bundles are endowed with 14,820 elements (14,821 DOFs) and 67,456 elements (67,457 DOFs), respectively.Given the space resolution of the biventricular mesh, we apply non-Gaussian quadrature rules to recover convergent conduction velocities [62].We consider a transmural variation of ionic conductances to differentiate epicardial, mid-myocardial and endocardial properties [63].To solve Eq. , we leverage an Implicit-Explicit time discretization scheme, where we first update the variables of the ionic model and then the transmembrane potential [13].Specifically, in the monodomain equation, the diffusion term is treated implicitly and the ionic term is treated explicitly.The latter is discretized by means of the Ionic Current Interpolation scheme [28].We prescribe the fiber distribution according to a Laplace-Dirichlet Rule-Based algorithm with α epi = −60

Branched Latent Neural Maps
We construct a geometry-specific surrogate model of cardiac function by building a feedforward partially-connected NN that explores the variability of our physics-based electrophysiology model detailed in Section 3.2 while structurally separating the role of temporal t and functional θ EP parameters.This recent scientific machine learning tool, proposed in [53], allows for different levels of disentanglement between inputs and outputs.The surrogate model reads: Weights and biases w ∈ R Nw encode the algebraic structures of a feedforward partially-connected NN, which represents a map BLNM : R 1+N P → R Nz from time t and N P = 7 cell-to-organ scale electrophysiology parameters 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), aV L(t), aV R(t) and aV F (t) a posteriori by means of Equations ( 3) and (4).Furthermore, vector z(t) leverages some z latent (t) latent variables that enhance the learned temporal dynamics by acting in regions with steep gradients [53].
We perform nonlinear optimization with the BFGS algorithm to tune NN parameters.In particular, we monitor the MSE of surrogate vs. physics-based ECG pseudopotentials to find an optimal set of weights and biases w, that is: where z leads (t) ∈ [−1, 1] 9 represents BLNM outputs and z numerical (t) ∈ [−1, 1] 9 defines the physicsbased numerical simulations, both in non-dimensional form.Time t ∈ [0, 1] and model parameters θ EP ∈ [−1, 1] N P are also normalized during the training and testing phases.We refer to [53] for a detailed description of all the properties related to BLNMs that enable them to effectively learn complex physical processes.

Parameter estimation
We employ our trained BLNM to find a set of model parameters θ EP that matches z ECG (t) ∈ R 12 with z clinical (t) ∈ R 12 .Here, z ECG (t) is the vector of BLNM physical outputs z leads (t) manipulated according to Equations ( 3) and ( 4) to generate the full 12-lead ECGs, and z clinical (t) is the clinically measured ECG data vector.In particular, we perform derivative-free optimization by employing the Nelder-Mead method [33], where we specify a loss function given by the MSE of the mismatch between the trained surrogate vs. clinical ECG potentials, that is: which leads to a set of calibrated model parameters θ NM EP and corresponding 12-lead ECGs z NM ECG (t).We initialize our optimization algorithm with a random set of model parameters θ init ∈ [−1, 1] N P .We repeat the optimization process 100 times and we average the model parameters obtained during the different trials in order to get θ NM EP .

Sensitivity analysis
We perform a variance-based sensitivity analysis using Shapley effects [56] 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 [12] to define the input multivariate distribution, which is given by a Gaussian copula and a series of N P marginals defined by standard normal distributions centered in θ NM EP , that is N θ NM,i EP , 0.2 for i = 1, ..., N 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 (7).We fix 2000 permutations, 500 bootstrapped samples, and 50 samples to estimate conditional variance for 3 times.

Uncertainty quantification
We employ a BLNM within HMC [5] 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 P( θ EP |x), given a certain prior probability distribution P( θ EP ) with respect to the model parameters in non-dimensional form θ EP ∈ [−1, 1] N 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 [21].This algorithm, which shares and enhances some of the features of sequential [26] and differential evolution [61] 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 [21].For further details about the mathematical derivation of HMC and its application to cardiac simulations we refer to [54].
We run 4 chains with 1,000 adaptation samples in the warm-up phase and 1,000 effective samples to estimate the posterior distribution, with a fixed 90% acceptance rate.For all model parameters, we consider prior distributions where θ NM EP ∈ [−1, 1] N 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 ι = 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 modeling error of reducedorder 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 multivariate normal distribution centered in the BLNM predictions z ECG (t) for the given patient-specific observations z clinical (t), which reads: σ meas = 0.1 is the a priori fixed standard deviation dictating the measurement error [15,69], whereas is the exponentiated quadratic kernel of a zero-mean Gaussian process GP(0, k( t, t ′ ; σ GP , l GP )) [46].Amplitude σ GP ∼ N (0.01, 1.0) and correlation length l GP ∼ N (0.01, 1.0) are additional hyperpriors tuned during HMC to quantify the surrogate modeling error, which may change according to the specific observation.Vectors t and t ′ represent discrete time points in the [0, 1] interval.A full covariance matrix in the multivariate normal distribution allows us to model the correlation among different leads.We evaluate convergence of the HMC chains by checking that the Gelman-Rubin diagnostic provides a value less than 1.1 for all the model parameters θ EP , l GP and σ GP [6,17,65].

Software and hardware
All electrophysiology simulations are performed at the Stanford Research Computing Center using svFSIplus [70], a C++ high-performance computing multiphysics and multiscale finite element solver for cardiac and cardiovascular modeling.This solver is part of the SimVascular software suite for patient-specific cardiovascular modeling [64].
We train the NNs by using BLNM.jl[22,45], 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[3] for parallel hyperparameter optimization by combining the Message Passing Interface (MPI) with Open Multi-Processing (OpenMP) on physical and virtual cores, respectively.
We perform sensitivity analysis and parameter estimation with uncertainty quantification using GlobalSensitivity.jl[11] and Turing.jl[16], 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.

Discussion
We present a complete computational pipeline to build digital twins of cardiac electrophysiology for congenital heart disease in pediatrics.This cohort of patients is understudied in cardiology [32,62], as multiphysics and multiscale numerical simulations are mostly focused on adults with certain sets of pathologies, such as dilated, ischemic and hypertrophic cardiomyopathy, arrhythmias or bundle branch block [36,38,52,58].
In this pipeline, we leverage biophysically detailed and anatomically accurate computational electrophysiology models, a recently proposed scientific machine learning tool for surrogate modeling, and robust Bayesian inference methods for personalized calibration of model parameters to match clinical 12-lead ECGs of an HLHS pediatric patient.We certify the impact and reliability of our estimation against clinical recordings by integrating fast and effective sensitivity analysis and uncertainty quantification.We run electrophysiology simulations with the estimated model parameters in order to investigate different scenarios of clinical interest in silico.We conclude that this pediatric patient presents activation and repolarization patterns similar to a left bundle branch block, where the interventricular dyssynchrony and the geometrical personalization of the Purkinje network play a minor role with respect to conductances and conductivities, even for QRS complex calibration.
Image processing allows us to get all the anatomy-specific features of this pediatric patient and our calibration of cell-to-organ level model parameters enables patient-specific electrophysiology simulations.Nevertheless, given the non-convexity of the optimization problem, it is important to stress that the final set of model parameters might not be unique and there could be other choices that lead to similar approximation errors against clinical recordings.Indeed, we notice that changing random seeds or trying different optimizers, such as second-order local BFGS or even global Adaptive Differential Evolution [68], may have an influence on the initial parameter estimation.These options are available within the BLNM.jllibrary.However, these effects are accounted for and mitigated by averaging many different trials and by running uncertainty quantification.
Performing ad-hoc sensitivity analysis for a specific parameter calibration provides individualized information, as these assessments may change on a patient to patient basis.Furthermore, we underline that sensitivity and practical identifiability (or trustworthiness) of model parameters are generally correlated.For instance, the maximum rapid delayed rectifier current conductance G Kr and the level of interventricular dyssynchrony t stim LV have the lowest relative impact on this 12-lead ECGs personalization (see Figure 4) and present the highest degree of uncertainty, that is a wider posterior distribution, among all physics-based model parameters (see Figure 5).
While the computational pipeline encompasses several rigorous steps, the physics-based model still requires high-performance computing and longer computational times compared to other approaches for digital twinning on adults [19,20], which rely on more phenomenologically-based models but do not include robust methods for sensitivity analysis and uncertainty quantification.However, the monodomain equation, coupled with the ten Tusscher-Panfilov ionic model, provides an accurate mathematical model, where relevant model parameters with a direct physiological interpretation can be properly tuned.Moreover, its higher computational costs could be mitigated in future work by novel numerical methods in the framework of matrix-free [1] and Isogeometric Analysis [7].
A limitation of the presented approach lies in the lack of experimental validation of the parameter calibration process.Indeed, mathematical modeling of congenital heart disease requires several assumptions due to the current lack of information in pediatric populations regarding fiber orientation, Purkinje structure, ionic current conductances, and conduction velocities.Future studies should incorporate these data as they become available.Nevertheless, our estimations are robust, account for uncertainty quantification and are widely contained within the range explored by the electrophysiology simulations (see Table 1 and 3, along with Figure 5).
In future developments, we aim to encode anatomical variability and different CHDs, such as Tetralogy of Fallot, transposition of great arteries or atrial and ventricular septal defects within BLNMs.In this manner, the computationally expensive offline phase dictated by accurate numerical simulations and the training of the NN can be performed only once before being applied to new patients.Robust parameter estimation and uncertainty quantification will be then feasible for those CHDs within minutes, compatible with the time frame required by the clinical practice.
t e x i t s h a 1 _ b a s e 6 4 = " Z 0 y R c E k I c 0 4 b A E k 4 O f q / 6 t 7 y 6 k U = " > A A A D 2 X i c b Z N b b 9 M w F M f d h M s I t w 4 e e b G o m B B C V T u J y w v S x I q G t G o q 0 t p N a r r K c d z W q u 0 E 2 w F V I Q 8 8 g B C v f D P e + B R 8 B U 6 y r O s S L M X 5 n 3 N + P s c 5 d o J Y c G M 7 n T 8 N x 7 1 2 / c b N r V v e 7 T t 3 7 9 1 v b j 8 Y m S j R l A 1 p J C J 9 G h D D B F d s a L k V 7 D T W j M h A s J N g u Z / H T z 4 x b X i k j u 0 q Z h N J 5 o r P O C U W X N P t x l 9 f s c 8 0 k p K o M P V j o o l 8 N 8 h K N f U l s Q s t U 3 B l X p 0 E L o h E a F Y S X q l v F 8 y S K n i w T / p Z e r B O l Z t V 5 I h s E k f 1 H I d 6 E w C r B p g r g K k C P a j Q W 8 e J 4 j W A b w L c R D U g 3 g T i R C + r h O 2 P j O X Q F L v G + q P s 7 E I X s c z z / I D N u U r Z x 6 Q 4 g 2 f Z h Y c I P l c s B L s 8 B v w G j 3 e K B j 7 H e Z P y + V D D 3 M t l j 8 P k 4 x 3 Y G Y i y + O T s 2 P M Z 7 O Y y W W 5 d F p s 2 W 5 1 2 p x i 4 L r q l a K F y D K b N 3 3 4 Y 0 U Q y Z a k g x o y 7 n d h O U q I t p 4 J B / s S w m N A l m b M x S E U k M 5 O 0 u J k Z f g K e E M 8 i D Y + y u P B u r k i J N P n t A T L v k q n G c u f / Y u P E z l 5 P U q 7 i x D J F z w v N E o F t h P N r j k O u G b V i B Y J Q z W G v m C 6 g q 9 T C z + B B E 7 r V T 6 6 L 0 W 6 7 + 7 L 9 4 s N u a + 9 t 2 Y 4 t 9 A g 9 R k 9 R F 7 1 C e + g 9 G q A h o s 7 I + e J 8 c 7 6 7 Y / e r + 8 P 9 e Y 4 6 j X L N Q 3 R l u L / + A e 2 v T B 4 = < / l a t e x i t > ✓ EP = [G CaL , G Na , G Kr , D ani , D iso , D purk , t stim LV ] T r L 2 e e b 5 e e a t k 1 w K i 8 3 m 7 y A 8 d / 7 C x U t r l 6 M r V 6 9 d v 1 F b v 7 l n s 8 J w 6 P J

V6Figure 1 :
Figure1: Sketch of the computational pipeline.We reconstruct the patient-specific geometry with HLHS from imaging.We generate a dataset of electrophysiology simulations encompassing cell-toorgan variability in model parameters.We train a BLNM that effectively reproduces 12-lead ECGs while covering model variability.We employ the BLNM for digital twinning.

Figure 3 :
Figure 3: Physics-based electrophysiological modeling.Spatio-temporal membrane action potential evolution for one electrophysiology simulation in the dataset performed on the HLHS pediatric patient.

Figure 4 :
Figure 4: Sensitivity analysis for the seven model parameters encoded in the BLNM via Shapley effects.

Figure 5 :
Figure 5: Inverse uncertainty quantification: parameter uncertainty.One-dimensional views of the posterior distribution.Different colors represent different HMC chains.

Figure 7 :
Figure 7: Running an in silico clinical trial.Activation (top) and repolarization (bottom) maps with personalized model parameters (left), left bundle branch block (center) and right bundle branch block (right).

Table 1 :
Parameter space for cardiac electrophysiology sampled via latin hypercube for the numerical simulations performed with the physics-based mathematical model.

Table 2 :
Branched Latent Neural Map hyperparameter tuning.Original hyperparameter ranges and optimized hyperparameter values for the final training stage.

Table 3 :
Parameter estimation.Calibration with the optimized BLNM for cell-to-organ level model parameters of the physics-based mathematical model.The MSE between the BLNMs predictions and the clinical recordings, in non-dimensional form, is 0.16.