Introduction

Digital Twins are playing an increasingly important role in healthcare to enable tailored therapy design for precision medicine. The importance of this emerging technology was attested to by recent regulatory attention and the publication of the ASME V&V40 Standard in 2018¹. The Standard provided a risk-based framework for evaluating their credibility for biomedical applications.

At the core of the cardiac Digital Twin are realistic mechanistic models that can coherently assimilate multi-modal and multi-scale data and can explain disease mechanisms and predict therapy outcomes^2–4. An important aspect of cardiac function is the coupling of electrophysiology and mechanics in the ventricles^5–7. To describe electromechanical function, Digital Twin models are built on mathematical descriptions of electromechanical coupling at cellular, tissue, and organ scales. The biophysical detail gives the model the advantage of high explainability and predictive power over statistical or machine learning approaches.

Electromechanical cardiac Digital Twins have had broad applicability in cardiology, such as in unravelling the relationship between the electrocardiogram (ECG) and mechanical deformation⁸, the ECG and ejection fraction in post-myocardial infarction⁹, demonstrate how mechanical deformation provides triggers and substrate for arrhythmia¹⁰ and affects re-entrant wave stability¹¹, and predicting therapy outcomes in hypertrophic cardiomyopathy¹² and heart failure¹³. Recent developments have enabled increasingly personalised cardiac Digital Twins through assimilation of omics, imaging, blood pressure, and electrocardiogram (ECG) data^6,14–16 and promises to deliver patient-specific therapy planning¹⁷ in the near future.

Despite significant developments, a lack of thorough reviews of model credibility considering both electrophysiological and mechanical properties from cell to organ scales hinders wider acceptance of cardiac Digital Twin technologies. Thus far, cardiac applications of the ASME V&V40 standards have focused on device optimisation using fluid-structure interaction models for left ventricular assistive devices¹⁸ and for artificial heart valves¹⁹, but this has yet to be applied to ventricular electromechanical simulations. Examples of rigorous model evaluation at cellular scale^20–22 highlight the importance of strictly separating calibration and validation data²⁰ and consideration of both electrophysiological and mechanical biomarkers²³. At the ventricular scale, most studies focused either on the ECG⁸ or pressure-volume biomarkers^{24–27 18 28} but not on both, and deformation and strain biomarkers are rarely considered. Simulation studies have investigated how model parameters relate to specific aspects of cardiac electromechanical function, such as the atrioventricular plane displacement²⁸, and how model deformations relate to ECG morphologies⁸, but not in a comprehensive manner. Large scale global sensitivity analyses studies of pressure and volume biomarkers have analysed effects from anatomical variations²⁵ and contractility and haemodynamic parameter variations²⁹, without any analysis of concurrent ECG biomarker effects. Furthermore, the strategies and data used for calibration and validation in these studies were typically not discussed^5,6,25.

Therefore, in this study we aim to propose a calibration and validation strategy based on a thorough review of the literature for model parameters and considering both ECG and mechanical biomarkers for evaluating ventricular electromechanical Digital Twins, following ASME V&V40 principles. We apply this to a multiscale, biventricular model with the aim to 1) reproduce healthy pressure-volume, ECG, deformation, and strain dynamics, and 2) investigate multi-scale mechanisms that underpin the relationship between ECG and pressure-volume biomarkers.

Materials and methods

Multi-scale human biventricular electromechanics

Human electrophysiological model

Simulation of the electrical excitation and relaxation through the ventricles was through the monodomain model with orthotropic diffusion along the fibre, sheet, and sheet-normal directions. Ionic current dynamics and calcium handling dynamics were modelled using the human ventricular electrophysiological cell ToR_ORd model ²⁰, with extensive validation on experimental action potential morphology and calcium dynamics, as well as responses to multi-channel drug block.

Purkinje-myocardial junctions were modelled using a fast endocardial activation layer with isotropic conduction velocity σ _endo on the left and right ventricular endocardial surfaces. Standard 12-lead electrode positions were manually mapped from the heart and torso geometry of a previous study to the current geometry, and the 12-lead ECG was simulated at these electrode positions using the pseudo-ECG method¹⁵.

Mechanical model

The mechanical behaviour was characterised by the balance of linear momentum in a total Lagrangian formulation, considering inertial contributions but ignoring volumetric forces. The reference configuration was the scaled down version of the biventricular mesh (diastasis) Ω ⊂ R³, at t =0. The constitutive law for the passive mechanical behaviour of the ventricular tissue was a nearly incompressible version of Holzapfel and Ogden which was orthotropic in mechanical behaviour⁷. The strain energy density function defining this hyper-elastic continuum is:

where J=det (F), _C=F^T_F, bar notation indicates that , and indicate fibre and sheet directions, respectively, in the reference configuration such that , and K_ct controls the compressibility of the material. The reader is directed elsewhere⁷ for the precise form of the passive contribution to the second Piola-Kirchhoff stress tensor.

Boundaries were defined as: epicardium Γ_epi, left ventricular endocardium Γ_LVendo, right ventricular endocardium Γ_RVendo, and the non-endocardial surfaces of the valvular plugs Γ_valve,

where n₀ was the outward normal defined on the whole material boundary ∂Ω₀, and K_epi (uv c_l) was the stiffness corresponding to an elastic spring boundary condition applied to the endocardium and was described as a step function of the longitudinal coordinate l with magnitude of K_epi at values of 0<uv c_l <0.85 (i.e. applied for 85% of the longitudinal dimension of the heart from the apex). This was a simplified version of the method²⁸, which uses an exponential decay formulation at the ‘edge’ of the pericardial constraint rather than a step function. P_LVendo and P_RVendo were the pressures applied to the left and right endocardial surfaces, respectively, and was prescribed by two independent piece-wise functions that describe the five-phases (not including initiation) of the cardiac cycle:

• Phase 0: Initiation. Endocardial pressures were linearly increased to reach diastasis values P_RVDS and P_LVDS.

• Phase 1: Active diastolic filling. Endocardial pressures were linearly increased over time t _EDP to end diastolic values EDP_LV and ED P_RV to mimic the effect of atrial contraction.

• Phase 2: Isovolumetric contraction. Electrical activation ensues and active contraction develops. During this, the ventricular pressure was controlled such that the volume was maintained to be constant via:

  

  

  Where C _pLV, C _pRV, C_VLV, C_VRV, were the inverse of penalty terms for volume difference and volume rates for the left and right ventricles, respectively. The volume rate term was used as stabilisation against spurious oscillations of the ventricular pressure, which might occur due to inertial effects.

• Phase 3: Ejection. When the ventricular pressure surpasses the aortic P_AO and pulmonary artery P_PA pressures, the ejection phase begins. To model the blood pressure of the systemic circulation system during the ejection phase, the following two element Windkessel model was used:

  

  Where C _AO and C_PA were the compliance of the aortic and pulmonary arteries, respectively, and R_AO and R_PA

  

  were the resistance of the aortic and pulmonary circuits, respectively. During the ejection phase, the ventricular pressures P_LVendo and P_RVendo were modelled as equal to the arterial pressures P_AO and P_PA respectively, disregarding negligible pressure gradients across the arterial valves.

• Phase 4: Isovolumetric relaxation. This phase was triggered when the ventricular flow reverses V^• _LVendo>0, and V^•_RVendo>0, and it follows the same formulation as the isovolumetric contraction Phase 2.

• Phase 5: Passive filling. This phase begins when the endocardial pressure drops below the left and right atrial pressures P_LA, P_RA. The pressure was prescribed so that the two ventricular volumes were returned to its diastasis value (DSV_LV, DSV_RV).

Where C _pLAV, C _pRAV, C_vLAV, C_vRAV were the inverse of penalty terms for volume different to DSVs and volume rates, for the left and right ventricles, respectively.

The two ventricular volumes V _LVendo and V _RVendo were calculated and updated throughout the cardiac cycle by using the divergence theorem and the assumption that the volume was constrained by a flat lid located on the aortic and pulmonary valvular planes, for right and left ventricles respectively, leading to:

Where e _AO and e_PA were defined so that e _AO·n₀=0 on the aortic valve plane, and e_PA · n₀=0 on the pulmonary valve plane.

Electromechanical coupling

Electrophysiological activation begins at Phase 2 of the cardiac cycle and drives the active tension generation. This is modelled at the cellular level by coupling of the Land active force generation²¹ with intracellular calcium dynamics from the human cellular electrophysiology ToRORd model²⁰, henceforth referred to as the human cellular electromechanical ToRORd-Land model. In order to couple the ToRORd to the Land model, the Hill coefficient of calcium cooperativity and tropomyosin rate constant of the Land model were re-fitted in a previous study²³ so that physiological active tension could be produced in response to the calcium transient dynamics of the ToRORd model, which were different to the calcium transient to which the original Land model was fitted.

In the Land model, active force generation is modelled using four cross-bridge cycling states (blocked, unbound, pre-power stroke and the force-generating state) and the transition rates between each were calibrated to human skinned myocyte data²¹. The ToRORd-Land coupling is bidirectional, where calcium from the ToRORd binds to troponin C and drives the transition between the blocked and unbound states of the crossbridge system, while the binding of calcium from the Land model acts as a buffer for the ToRORd calcium transient. Furthermore, the length-dependence of troponin calcium affinity/sensitivity in the Land model means that calcium buffering in the ToRORd is dynamically affected by tissue deformation. The Land model was first constructed using skinned myocyte data, then several of the parameters were re-fitted to enable sufficient ejection at the ventricular scale, including the transition rate from the pre-power stroke to force-generating state of the cross bridge (k_ws), the sensitivity of troponin to calcium binding (Cal50)²¹. At the tissue level, the active tension was then added to the stress tensor to drive mechanical contraction⁷.

Literature-informed parameter initialisation

A review of the literature produced a list of initial model parameter values as compiled in Table 1, including also reported parameter variability. Passive mechanical properties were extracted from the non-viscoelastic version of the model fitted to human ex vivo stress-strain measurements³⁰. Ionic time constants and conductance values were set to baseline values of the human membrane kinetics ToR_ORd model, which had been calibrated to patch clamp, action potential and calcium transient data²⁰. The active tension Land model parameters were as previous calibrated to achieve a physiological calcium transient²³.

Table 1

Model assessment strategy

Following examples of the ASME V&V 40-based analysis applied in the context of cardiovascular science ^18,48, we devised a calibration and validation strategy for human ventricular electromechanical modelling and simulation frameworks by considering the following:

Question of interest

Can the modelling and simulation framework produce healthy physiological electromechanical function at cellular, tissue, and organ levels including both ECG and pressure-volume loop biomarkers?

Context of use

Human ventricular electromechanics models provide mechanistic explanations for various disease conditions and evaluate therapeutic interventions. Here we focused on model evaluation under healthy conditions to build the foundations of model credibility for disease and therapy investigations in the future.

Quantities of interest

To ensure model relevance to its context of use, the quantities of interest were chosen to be biomarkers with known links to key cardiac diseases, including Torsade de Pointes, myocardial infarction, hypertrophic cardiomyopathy, dilated cardiomyopathy, heart failure, and pulmonary hypertension. Therefore we selected the following biomarkers from the 12-lead ECG (QRS duration, QT dispersion, T wave amplitude, T peak to T end duration, and JT interval), the pressure-volume relationship (left and right end diastolic and end systolic volumes and end systolic pressures, and left ejection fraction, left pressure upstroke velocity, left peak ejection and peak filling flow rate), and from measurements of displacement (end diastolic and end systolic wall thicknesses, atrioventricular plane displacements, systolic percentage volume change, and apical displacement) and strain (systolic global longitudinal strain, radial strain, circumferential strain, and fibre stretch ratio). The relevance of each biomarker to cardiac diseases was summarised in Table 3.

Cellular biomarkers including action potential duration (APD90), active tension peak (Ta_max), active tension duration (TaD50) and peak upstroke of active tension (dTa/dt_max) were not explicitly included in the biomarker list since they have already been used to calibrate and validation the cellular electromechanical model previously^20,23. However, these properties were evaluated in uncertainty quantification (see below) to inform ventricular simulations.

Compilation of biomarker data

The healthy ranges for each quantity of interest (biomarker) were compiled from the literature. Reports from larger datasets were preferred for better estimates of population variability. For this reason, values were preferentially included from the UK Biobank, which is a large-scale biomedical database with health information from half a million UK participants, with more than 40,000 cardiac magnetic resonance recordings available^49–51. This is then supplemented by smaller studies as needed. The use of gold standard magnetic resonance imaging data was preferentially included over other imaging modalities such as echocardiography for volumetric and deformation biomarkers.

Separation of calibration vs validation data

The quantities of interest were separated into a calibration and a validation dataset²⁰. The calibration dataset consisted of ECG and pressure-volume biomarkers that clinically defined healthy electromechanical function, while the validation dataset included displacement and strain biomarkers that explored the model’s capabilities for providing mechanistic explanations.

Sensitivity analyses strategy

The model was further validated through sensitivity analyses. To showcase versatility, a wider set of parameters were included in the analyses than has explicitly been linked to any disease or therapeutic target.

The following model parameters were included from the cellular electrophysiology (conductances of the fast sodium current GNa, late sodium current GNaL, L-type calcium current GCaL, transient outward potassium current Gto, rapid delayed rectifier potassium current GKr, slow delayed rectifier potassium current GKs, sodium calcium exchanger GNCX, inward rectifier current GK1, sodium-potassium pump current GNaK, peak calcium release Jrel, peak calcium re-uptake SERCA, and magnitude of the stimulus current Istim), the cellular active tension generation (peak active tension T_ref, calcium sensitivity Cal50, cross bridge transition rates from unbound to pre-powerstroke k_uw and from pre- to post-powerstroke k_ws), passive mechanical properties (bulk stiffness a, fibre stiffness a_f, sheet stiffness a_s, fibre-sheet shear stiffness a_fs, compressibility coefficient Kct, pericardial constraint stiffness K_epi), and circulatory haemodynamics (aortic resistance R_LV, aortic compliance C_LV, aortic ejection threshold pressure P_ejectionLV).

Model parameter ranges were informed by variabilities in literature reports (summarised in Table 2). Cellular conductances were uniformly varied from 50 to 200% to cover the experimentally measured ranges as done in previous studies²². For passive mechanics, only the parameters carrying units of stiffness were included in the exploration and their ranges come directly from the literature. The compressibility parameter (Kct), pericardial constraint stiffness (Kepi), and the active tension coefficient (Tref) parameters were allowed greater ranges of variation since they do not correspond to experimental measurements and their physical meanings were model-dependent. Similarly, the resistance and compliance parameters and the ejection threshold pressure were allowed a large uncertainty range since they were part of a lumped parameter model and do not directly map to experimental measurements.

Table 2

A specific example of model calibration, validation and sensitivity analysis

Input data

To demonstrate our strategy we selected a patient-specific (female, age 49, weight 90 kg) biventricular tetrahedral mesh with valvular plugs from four-chamber meshes available from a published database²⁵. Since the mesh was constructed based on imaging data at end diastole, the mesh was scaled down isotropically to achieve a diastasis volume of 80 mL³⁹ to mimic the process of unloading to diastasis. Fibre f ₀, sheet s₀, and sheet-normal n₀ vectors at resting were generated using a rule-based method designed for biventricular geometry with outflow tracts⁴². Ventricular coordinates were taken from the published database²⁵, and consisted of the longitudinal coordinate (uv c_l) where uv c_t=1 at the apex and 0 at the base, a transmural coordinate (uv c_t), where endocardium was uv c_t=1 and epicardium was uv c_t=0, and a rotational coordinate (uv c_r), where uv c_r=0 was in the middle of the left ventricular lateral wall.

Calibration to ECG and pressure-volume biomarkers

The 12-lead ECG signals of a female healthy volunteer with similar myocardial mass (108 g) as the biventricular mesh (96 g) was extracted and beat-averaged¹⁵ for calibration of electrophysiological activation and repolarisation characteristics. The conduction velocities along the fibre and sheet-normal directions were extracted from surgical measurements⁵². The root nodes of earliest activation on the endocardial surface as well as the conduction velocities on the endocardial surface and transmurally were calibrated using the QRS segment of the ECG signal, while spatial heterogeneity in the slow rectifier potassium current (IKs) was calibrated to the ST-T segment of the ECG signal, using a Bayesian-based inference method^15,53. The calibrated model had a transmural conduction velocity of 48 cm/s, and the fibre and sheet-normal conduction velocities were set to 67 cm/s and 44 cm/s as measured using plunge needle in control subjects⁵².

When the initialised parameter values were used in a preliminary electromechanical simulation, the results undershot healthy values of left ventricular peak systolic pressure by 2 kPa, ejection fraction by 32%, stroke volume by 29 mL, and peak ejection rate by 142 mL/ms, and peak filling rate by 207 mL/ms, and overshot the end systolic volume by 18 mL, and the peak pressure upstroke velocity by 70 kPa/ms. This demonstrated the need for model calibration and explorations of model uncertainties.

Due to the high computational cost of ventricular electromechanics simulations, calibration methods that require high volume of evaluations were infeasible. In this study, we designed a sequential calibration strategy based on knowledge of model sensitivities (described in detail in Appendix 1). This involves a series of sampling and parameter selections as follows, where single parameters were uniformly sampled and multiple parameters were sampled using Latin Hypercube Sampling:

1. Sample ejection pressure threshold, Kct, Cal50 and k_ws and select parameter set that gives highest LVEF.

2. Sample arterial resistance R_LV and select value that gives best peak systolic pressure.

3. Sample GCaL and select value that gives highest LVEF.

4. Sample k_ws and select value that gives best peak ejection rate, dP/dt_max and LVEF.

5. Sample diastolic volume change parameter and select value that gives best peak filling rate.

Validation by comparison to deformation and strain biomarkers

The calibrated model was evaluated against displacement and strain biomarkers as listed in Table 3 and not used for model calibration. Atrioventricular plane displacement was calculated by tracking the longitudinal displacement of the top 10% of the geometry using the apex-to-base ventricular coordinate. Wall thickness was calculated using perpendicular projections of the endocardial and epicardial surfaces. Myocardial strains were evaluated at a mid-ventricular short axis slice and at a four-chamber longitudinal, to match DENSE measurements.

Circumferential, radial, and fibre strains were evaluated using the short axis slice, while the longitudinal strain was evaluated using the long axis slice. The strain transients were then time-shifted to begin at the end of diastolic filling (Phase 1) and offset by the end diastolic strain, to match DENSE measurements.

Sensitivity analysis at cellular and ventricular scales

A global sensitivity analysis was first performed at the cellular scale using the human cellular electromechanical ToR_ORd-Land model²³ to identify the key parameters that affect active tension and action potential duration. The conductances of all ionic currents as well as troponin calcium sensitivity and the cross-bridge power stroke rate and myosin head attachment rates (k_ws and k_uw) were varied from 50% to 200%. The key parameters that affect active tension dynamics and action potential duration were identified through ranking the parameters using the total Sobol index. This analysis showed that the key parameters which affected active tension amplitude, duration, and upstroke were GCaL, SERCA, calcium sensitivity Cal50, the cross-bridge cycling rate (k_ws) and GNaL (Appendix Figure 1).

Figure 1

At the ventricular scale, the mechanical and haemodynamic parameters listed in Table 2 along with the cellular parameters GCaL, SERCA, Cal50, kws, and GNaL were included one-at-a-time sensitivity analyses, varying each parameter uniformly in a total of eight simulations each over the ranges as specified in Table 2 for mechanical and haemodynamic parameters and over 50% and 200% for cellular parameters. Each biomarker as listed in Table 3 were evaluated, and linear regression was performed.

Parameter-biomarker relationships that had p-values less than 0.05 and absolute r-values larger than 0.6 were considered significant and linear relationships. The magnitudes of biomarker effects were normalised by the maximum magnitude for each biomarker. To visualise the relative importance of model parameter on simulated biomarker, a summary diagram was generated with lines linking each parameter to each biomarker. The colour of each line indicated a positive (red) or a negative (blue) relationship. The thickness of each line was scaled by the normalised magnitude of the effect. The transparency of each line was scaled by the r-value as measure of the linearity of relationship (Results Figure 3).

Computations and software

Cellular electrophysiological simulations and Latin Hypercube Sampling were performed using bespoke MATLAB codes. Coupled cellular electromechanics, as well as biventricular electromechanics simulations, were performed using the high-performance numerical software, Alya for complex coupled multi-physics and multi-scale problems⁵⁴ on the JURECA pre-exascale module supercomputer operated by Juelich Supercomputing Centre, Germany, through a PRACE-ICEI (project, as well as on the ARCHER2 supercomputer provided by the UK National Supercomputing Service through the CompBioMedX project (Computational Biomedicine a the Exascale) funded by EPSRC under grant agreement EP/X019446/1. Code for the generation of simulation files for multi-scale sensitivity analysis, calibration and validation, as well as scripts for evaluation of ECG, PV, deformation and strain biomarkers from ventricular simulations can be found at: https://github.com/jennyhelyanwe/Alya_input_setup/

Results

Experimental and clinical dataset for model calibration and validation

Table 3 presents the experimental and clinical dataset for healthy human ventricular electromechanics. For pressure-volume biomarkers and ECG biomarkers, reported values from the UK Biobank¹ were used, where the number of samples reaches the tens of thousands. Catheter measurements of intraventricular pressure dynamics were rare in the literature and have small dataset sizes. Mechanical displacement and strain biomarkers came from MRI studies with specialised sequences (e.g. DENSE strain imaging).

Table 3

Biventricular electromechanics model calibration to ECG and pressure-volume biomarkers

Following model calibration, the simulated ECG showed good R wave progression from V1 to V6 and QRS duration that fell within the ranges shown in Table 1 (Figure 1A) as well as simulated left ventricular ejection fraction of 51% and peak systolic pressure of 14 kPa, with end diastolic and end systolic volumes of 110 mL and 50 mL (Figure 1B), and a peak filling rate of 310 mL/ms (Figure 1C), all falling within 95% interval of physiological biomarker ranges as listed in Table 1. The right ventricular end diastolic and end systolic volumes also fell within physiological ranges, though the right ventricular ejection fraction did not, as this was not a key target for calibration. The calibrated model overshot the peak pressure upstroke by 30 kPa/ms and peak ejection rate by 180 mL/ms (Figure 1D). However, the importance of matching to these two biomarkers is lower compared with other haemodynamic markers due to the relative paucity of reported data: n=17 for peak pressure upstroke and n=15 for peak ejection rate, compared with n=804 for ejection fraction and n=150 for peak systolic pressure (Table 3), It is also possible that the calibration could be improved by including some additional model parameters, such as those relating to stretch-rate dependent force generation, even though they did not show a strong effect on single cell active tension upstroke (Appendix Figure A1). These additional parameters might have a stronger effect on the upstroke at ventricular scale than at single cell scale because the single cell simulations were performed under isometric conditions. The calibrated model parameters are listed in Figure 1E.

Validation by comparison to deformation and strain biomarkers, not used in calibration

Simulation showed broad agreement with the trend and magnitude of in vivo strain measurements from DENSE+cDTI images⁴⁷ for fibre strain, circumferential, radial, and longitudinal strain in mid-ventricular and four-chamber view slices. Simulated strain traces (n=5004 for short axis, n=9282 for long axis) showed higher variability than in vivo traces (n=30). Simulated fibre strain (Eff, Figure 2A) showed agreement in terms of transmural homogeneity when compared with in vivo measurements, and the magnitude of fibre shortening was similar to in vivo (compare with red line at −0.15 strain). Peak circumferential (Ecc), radial (Err), and longitudinal strains (Ell) in the simulations (Figure 2B) showed good agreement in magnitude with in vivo measurements (compare with red line at ±0.15 strain). The timing of simulated peak strains occurred earlier than that in vivo. However, this is partially because the model was not calibrated to the cardiac cycle timings of the single individual in the DENSE+cDTI dataset.

Figure 2

Figure 3

In terms of deformations, the simulated systolic vs diastolic atrioventricular plane displacement of the calibrated biventricular (1.5 cm) showed good agreement with ranges reported in the literature (1.4 to 1.9 cm, Table 1).

When compared with literature reports of wall thickness changes, simulations showed thickening in systole in agreement with literature. The simulated end diastolic wall thickness (0.5 cm) was lower than the report range [0.57, 1.38] cm (Table 1, n=4620), and the end systolic wall thickness (0.7 cm) was also lower than the report range of [1.07, 1.75] cm (Table 1, 4620). The low wall thickness in both diastole and systole were consistent with the small heart size of the example used (96 g ventricular mass), and a larger and thicker heart would have greater wall thickening and fall better within the population ranges.

One-at-a-time sensitivity analysis

Figure 3 summarises how variability in single model parameters relate to simulated ECG (Figure 3A), pressure-volume (Figure 3B), displacement (Figure 3C), and strain (Figure 3D) biomarkers. Positive relationships between parameter and biomarker were shown in red and the negative were shown in blue.

ECG biomarkers (Figure 3A) were mostly affected by variabilities in ionic conductances (red block), with some small effects from the mechanical parameters (green block). QRS effects were very minimal, since conduction velocities were not altered in the analysis. QRS duration altered by SERCA (magnitude of change: 2.8 ms), GCaL (1 ms), GKr (1.8 ms), and GNaL (1.3 ms), while normalised QRS amplitude in lead V3 was very weakly affected by calcium sensitivity Cal50 (magnitude of change: 0.04), pericardial stiffness k_epi (0.01), and the passive stiffness parameter a_f (0.02) (Appendix Figure A1.A). The ST and T wave segments were much more sensitive to model parameters. The QT interval was strongly affected by GKr (magnitude of change: 173 ms), GNaL (92 ms), GNaK (72 ms), GCaL (69 ms), and SERCA (49 ms), with weak effects from calcium sensitivity (Cal50) (1.7 ms).

Normalised T wave amplitude was affected by GKr (magnitude of change: 0.1), GNaK (0.06), GNaL (0.06), and SERCA (0.04), GCaL (0.03), with very minor (<0.01) amplitude changes in select precordial leads for Tref (lead V4), Cal50 (lead V4), pericardial stiffness (lead V3) and Kct (lead V2) (Appendix Figure A1.B).

Pressure volume biomarkers (Figure 3B) were mostly affected by the mechanical (green block) and circulatory parameters (yellow block), alongside ionic currents (red block) that affect the calcium signaling system (GCaL and SERCA). Appendix Figure A2 illustrates in detail the parameter effects on the pressure-volume loops. The ejection fraction was predominantly affected by cross-bridge power-stroke rate (k_ws) (magnitude of change in ejection fraction: 32 %), compressibility (Kct) (16 %), GCaL (15 %), peak active tension (Tref) (12 %), SERCA activity (10 %). There were also some effects on LVEF from GNaL (9%), arterial resistance (6 %), calcium sensitivity (Cal50) (8 %), and very weak effects from pericardial stiffness (4%), passive mechanical parameters a (3 %), a_f (2 %), a_s (4 %), and GNaK (2%), and GKr (2 %).

The peak systolic pressure was affected by the same parameters as the ejection fraction, except for GNaK, GNaL and GKr, and with the addition of arterial compliance (magnitude change: 1.1 kPa) and ejection pressure threshold (1.9 kPa). The peak rate of rise of ventricular pressure (dPdtmax) was strongly influenced by only the cross-bridge power-stroke rate (k_ws) (magnitude of change: 1 kPa/s), with weak effects from calcium sensitivity (Cal50) (0.6 kPa/s) (Appendix Figure A3.A). The peak ejection flow rate (PER) was influenced by Cal50, Tref, k_ws, while the peak filling rate (PFR) was influenced by a larger group of parameters including Tref, Cal50, Kct, GCaL, SERCA, and the bulk parameter controlling the diastolic filling rate (Appendix Figure A3.B).

Mechanical displacement biomarkers (Figure 3C) were predominantly affected by active mechanics (green block) and calcium system conductances (red block), with very minor effects from passive mechanical parameters. Systolic AVPD (Appendix Figure A4.A) was strongly influenced by Cal50 (magnitude of effect: 0.54 cm), GCaL (0.36 cm), Tref (0.33 cm), Kct (0.41 cm), SERCA (0.31 cm), with weaker effect from k_ws (0.15 cm), arterial resistance (0.09 cm), and pericardial stiffness (0.08 cm). Systolic wall thickness was predominantly affected by the same parameters as systolic AVPD, except for Kct, which had only a weak effect. Systolic percentage volume change (Appendix Figure A4.C) was strongly affected by the compressibility parameter (Kct) (magnitude of change: 12 %), as was expected, but also strongly affected by GCaL (9 %), Tref (8 %), SERCA (6 %), and kws (3 %). An increase in SERCA and Cal50 shifted the timing of peak volume reduction to earlier during systole.

Simulated strain biomarkers (Figure 3D) were predominantly affected by passive mechanical parameters (green block) and circulatory parameters (yellow block). The influence of parameters on the radial, circumferential and longitudinal strain patterns (Appendix Figure A5) followed broadly that of the AVPD and wall thickness. Of note was that at high values of the ejection pressure, the mid-ventricular circumference expands rather than contracts during systole and fails to return to the reference condition even after relaxation, pointing to failure of proper pumping function (Appendix Figure A5.A).

Multi-scale mechanisms underpinning ECG and LVEF relationship

Simulation results indicated that the link between ECG biomarkers and LVEF was chiefly underpinned by GCaL and SERCA, which had strong effects on both T wave biomarkers and the pressure-volume loop (Figure 3). Mechanical parameters T_ref, calcium sensitivity (Cal50), compressibility (Kct), and fibre stiffness (a_f) had much stronger effects on LVEF than on ECG signatures. On the other hand, ionic conductance parameters GNaL, GKr, GNaK had much stronger effects on the T wave than on LVEF.

Figure 4 illustrates the multi-scale mechanisms of GCaL’s effect on both LVEF and ECG. A four-fold increase in GCaL caused an eight-fold increase in cellular active tension peak and a 200 ms prolongation of active tension duration (Figure 4A), which resulted in a 15 mL decrease in end systolic volume and 2.5 kPa higher peak systolic pressure at the ventricular scale (Figure 4C).

Figure 4

GCaL’s dramatic effect on end systolic volume was facilitated by 10 % higher systolic myocardial volume compression, 4.5 % higher systolic longitudinal contraction, 3.5 % higher systolic circumferential contraction, and 0.4 cm lower (towards apex) systolic position of the atrioventricular plane (Figure 4B). The increase in GCaL caused an increase of 60 ms in cellular action potential duration (Figure 4A), which prolonged the global repolarisation time (Figure 4D) and led to the prolongation of the QT interval (Figure 4C). GCaL did not have a significant effect on the activation pattern and conduction velocity in our simulations (Figure 4D) or on the QRS complex of the ECG (not shown).

Figure 5 illustrates how a four-fold increase in SERCA activity caused a 12% decrease in LVEF, with non-monotonic effects on the QT interval. The increase in SERCA conductance increased diastolic residual active tension (Figure 5A), which was sufficient to cause 3 mL reduction in end diastolic volume (Figure 5C). Even though increasing SERCA increased peak active tension at the cellular scale by 9 kPa (Figure 5A), it reduced peak systolic pressure and LVEF at the ventricular scale (Figure 5C). This was because increasing SERCA hastened the arrival of peak active tension by 50 ms and reduced its duration by 100 ms (Figure 5A), which meant that the mechanical contractions were not given enough time to developed, resulting in earlier and weaker systolic contractions (Figure 5B). On the electrophysiological side, while SERCA did not elicit strong APD changes at the single cell simulation, when combined with non-isometric conditions in ventricular simulations, it caused action potential prolongation at values higher than 1.5-fold and resulted in prolonged repolarization (Figure 5E).

Figure 5

Discussion

In this study we presented the calibration, validation and sensitivity analysis of human healthy ventricular electromechanical modelling and simulation, following the ASME V&V40-based assessment guidelines. The first step was compiling biomarkers values from multi-modal data characterizing pressure-volume, ECG, displacement, and strain behavior of human healthy ventricles, and separated the data for calibration and validation. The calibration was informed by unravelling the interlinking relationships between simulated biomarkers and model parameters, especially highlighting the plethora of model parameters that affect the ejection fraction. We demonstrated multi-scale mechanistic explanations of the relationship between LVEF and ECG underpinned by variability in L-type calcium channel conductance and SERCA activity.

Previous examples of applying the ASME V&V40 framework for cardiac model credibility assessment have focused primarily on evaluating computational models of medical devices such as the left ventricular assist device¹⁸ and artificial heart valves¹⁹, where model validation was performed by comparing to experimental recordings specifically designed for this purpose. These studies, as well as other examples of model validation without explicit reference to the ASME Standard, commonly have narrow focuses on specific diseases and/or therapies, and the biomarkers used for model evaluation were similarly limited in scope. In the case of disease models, calibrations were commonly performed at the cellular or tissue scale and validation was performed at the ventricular scale through comparisons with known ECG or ejection fraction phenotypes^9,12. While these studies have been valuable for providing model credibility in their specific disease/therapeutic contexts of use, they have not been designed to give credibility to the underlying computational electromechanics framework. In this study, we instead focused on building a comprehensive calibration and validation strategy for baseline healthy biventricular electromechanics. Through considering biomarkers at multiple scales and across both electrophysiology and mechanics, we hoped to establish a firm foundation of credibility for computational electromechanics that can be refined in future studies.

We showed in this study that literature values for model parameters do not automatically provide good baseline simulations. We proposed a sequential calibration strategy based on knowledge of model sensitivities, as detailed in the methods. This approach can be adapted for more sophisticated numerical techniques involving the use of emulators and Bayesian statistics to enable faster and more comprehensive explorations of the parameter space^24,26,86.

Our sensitivity analysis results showed that in the simulated ECGs, the T wave was far more sensitive to electrophysiological parameters than to mechanical ones, within the ranges of uncertainty expected from literature reports of each parameter. In our simulations, the ionic conductances GNaK and GNaL primarily affected the T wave through altering the action potential duration (Appendix Figure A1) while having minimal effect on mechanical deformation or strain. GCaL and SERCA strongly affected the repolarization pattern and deformation and strain patterns during systole (Figures 4 and 5). However, other parameters that had a strong effect on deformation and strain, such as Tref, Cal50, pericardial stiffness and Kct had only a minor effect on the T wave (Appendix Figure A2). Previous studies on this question focused on comparing either end diastolic versus end systolic geometry⁸ or a static versus dynamic geometry⁸⁷ when simulating the T wave. These simulations involved a much larger perturbation in deformation than relevant to physiological situations and showed a larger effect on the T wave amplitude than our simulations. Within the context of physiological systolic deformations, however, our results showed that deformation uncertainties had only minor effects on T wave amplitude.

Our simulations showed only minor changes to QRS amplitudes in response to variation in those parameters that influenced the end diastolic geometry: calcium sensitivity, through altering diastolic residual active tension, pericardial stiffness, through restricting diastolic inflation, and the passive mechanical parameters a and af, through altering bulk stiffness and stiffness along the myofibre direction. However, it is possible that the magnitude of changes to the QRS was underestimated in our simulations, because the local activation times on the endocardial surface in our simulations were fixed. This meant that our simulations did not explore the effect of variations in end diastolic geometry on the conduction velocity of the fast endocardial activation layer or on Purkinje conductivity. It should also be noted that other parameters that affect the activation pattern, such as GNa and conduction velocities, were not included in this analysis.

The fact that the initialised literature values of model parameters failed to achieve LVEF above 50 % pointed to limitations in using ex vivo tissue measurements to represent in vivo function. In our simulations, we saw that the LVEF was strongly sensitive to changes in the incompressibility of the tissue (Kct), such that an increase in compressibility of the myocardial tissue helped to increase LVEF. The systolic volume has been reported to change up to 13%⁸¹, pointing to limitations of assumptions of incompressibility in the modelling of the myocardial tissue. Additionally, our model did not achieve a high amount of wall thickening in validation, which meant that the high compressibility required in our model to achieve LVEF could be compensating for a lack of sufficient wall thickening, in addition to anatomical effects. Previous experimental and simulations studies have shown how wall thickening is related to sliding of sheetlets^88,89. However, in our simulations we did not see a significant effect of the sheer stiffness parameter a_fs on wall thickness or LVEF. One possible avenue of exploration would be the effect of sheet and fibre orientations on the thickening effect, as explored in a previous simulation study⁸⁹.

Displacement and strain biomarkers (Figure 2D) were far more sensitive to calcium handling dynamics, cross-bridge cycling rate, and contractility than to passive stiffness parameters, which indicates they would be better indicators of systolic rather than diastolic dysfunctions in HF. The sensitivity of strains and displacements to the pericardial stiffness parameter was unsurprising due to the association of reduced strains in diseases affecting the pericardial sack⁹⁰. The sensitivity of longitudinal and radial strain to the incompressibility of the myocardium (Kct) mirrors the sensitivity of the LVEF to that parameter, which was in accordance with the fact that those two strains were surrogate markers of LVEF and show better predictive power^83,84.

Conclusions

In this study, we show the importance of applying vigorous VVUQ evaluations of ventricular electromechanics for achieving physiological simulations. We set the basis for a strategy for calibrating and validating baseline electromechanical modelling and simulation frameworks for cardiac Digital Twins. We compiled a comprehensive list of biomarkers for evaluating healthy electromechanical function, and grouped the dataset into non-overlapping calibration and validation sets. We provide evidence of credibility through calibration to haemodynamic and ECG biomarkers and validation by comparison to strain and deformation biomarkers and uncertainty quantification. Furthermore, our analyses highlighted the interplay between cellular, tissue, and haemodynamic parameters on the LVEF and provided multi-scale explanations of its link with ECG biomarkers. Taken together, the study paves the way towards credible electromechanical cardiac Digital Twins.

Appendix

1. Pressure volume calibration strategy design

The goal of the calibration of the initialised set of parameters was to increase LVEF, increase systolic pressure, increase peak filling rate, but decrease peak ejection rate and lower dP/dt_max. The LVEF and peak systolic pressure had the highest degree of importance for calibration, since they were implicated in a variety of diseases and was well established. The filling rates and dPdt_max took a secondary role, since the data came from much smaller sample sizes, and the filling rates were measured using echocardiography techniques, which were not as reliable as volume measurements in CMR (Table 1).

From the uncertainty quantification, we saw that a handful of parameters can increase both LVEF and peak systolic pressure simultaneously: Tref, k_ws, GCaL, SERCA, Cal50, Kct (Figure 3). However, Tref and kws effects saturate at high values, large changes in SERCA and GCaL can become arrhythmic substrates, Cal50 has a non-monotonic effect on LVEF due to its effect on the diastolic function, and dramatic reductions in Kct brings unrealistic amounts of systolic volume change. A combination of changes in these parameters was needed to dramatically improve LVEF and peak systolic pressure, and in practice, was insufficient by themselves. Decreasing arterial resistance was also needed to help increase LVEF. However, this comes at the cost of reducing peak systolic pressure, which needed to be counterbalanced by increasing the ejection pressure threshold. An increase in ejection pressure increased time spent in isovolumic contraction and isometric force development, and thereby increased the peak systolic pressure. An additional challenge was that there was not a single parameter in our analysis that could decrease peak ejection flow rate (Tref, kws, Kct) and dPdt_max (kws) without also decreasing LVEF at the same time. This meant that we needed to achieve higher than healthy LVEF first. With all this in mind, the calibration strategy was as detailed in the Methods section.

2. Cellular global sensitivity analysis

Figure A1:

3. Additional ventricular sensitivity analysis results

Figure A2

Figure A3:

Figure A4:

Figure A5:

Figure A6:

Acknowledgements

This work was funded in whole, or in part, by the Wellcome Trust (214290/Z/18/Z). For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. This work was supported by a Wellcome Trust Fellowship in Basic Biomedical Sciences to B.R. (214290/Z/18/Z), the Personalised In-Silico Cardiology (PIC) project, the CompBioMed 1 and 2 Centre of Excellence in Computational Biomedicine (European Commission Horizon 2020 research and innovation programme, grant agreements No. 675451 and No. 823712), the Oxford BHF Centre of Research Excellence (RE/13/1/30181), PRACE-ICEI funding projects icp005, icp013, icp019.