Next Article in Journal
Fractional Growth Model with Delay for Recurrent Outbreaks Applied to COVID-19 Data
Next Article in Special Issue
An Automata-Based Cardiac Electrophysiology Simulator to Assess Arrhythmia Inducibility
Previous Article in Journal
A Novel βSA Ensemble Model for Forecasting the Number of Confirmed COVID-19 Cases in the US
Previous Article in Special Issue
Sensitivity Analysis of In Silico Fluid Simulations to Predict Thrombus Formation after Left Atrial Appendage Occlusion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Integrated Workflow for Building Digital Twins of Cardiac Electromechanics—A Multi-Fidelity Approach for Personalising Active Mechanics

1
Gottfried Schatz Research Center for Cell Signaling, Metabolism and Aging—Division of Biophysics, Medical University Graz, 8010 Graz, Austria
2
NAWI Graz, Institute of Mathematics and Scientific Computing, University of Graz, 8010 Graz, Austria
3
BioTechMed-Graz, 8010 Graz, Austria
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(5), 823; https://doi.org/10.3390/math10050823
Submission received: 31 January 2022 / Revised: 24 February 2022 / Accepted: 28 February 2022 / Published: 4 March 2022

Abstract

:
Personalised computer models of cardiac function, referred to as cardiac digital twins, are envisioned to play an important role in clinical precision therapies of cardiovascular diseases. A major obstacle hampering clinical translation involves the significant computational costs involved in the personalisation of biophysically detailed mechanistic models that require the identification of high-dimensional parameter vectors. An important aspect to identify in electromechanics (EM) models are active mechanics parameters that govern cardiac contraction and relaxation. In this study, we present a novel, fully automated, and efficient approach for personalising biophysically detailed active mechanics models using a two-step multi-fidelity solution. In the first step, active mechanical behaviour in a given 3D EM model is represented by a purely phenomenological, low-fidelity model, which is personalised at the organ scale by calibration to clinical cavity pressure data. Then, in the second step, median traces of nodal cellular active stress, intracellular calcium concentration, and fibre stretch are generated and utilised to personalise the desired high-fidelity model at the cellular scale using a 0D model of cardiac EM. Our novel approach was tested on a cohort of seven human left ventricular (LV) EM models, created from patients treated for aortic coarctation (CoA). Goodness of fit, computational cost, and robustness of the algorithm against uncertainty in the clinical data and variations of initial guesses were evaluated. We demonstrate that our multi-fidelity approach facilitates the personalisation of a biophysically detailed active stress model within only a few (2 to 4) expensive 3D organ-scale simulations—a computational effort compatible with clinical model applications.

1. Introduction

Cardiovascular diseases accounted for 32% of global deaths in 2019 and remain the leading cause of death worldwide [1]. Improvements in diagnosis, therapy stratification, and planning are needed to offer personalised precision therapies that are tailored to individual patients. Computer modelling of cardiac functions show promise in this regard, owing to its unique ability of integrating disparate clinical data into a quantitative and mechanistic framework that facilitates, ultimately testing, the prediction of outcomes for various therapeutic options [2,3,4,5].
Such advanced modelling applications rely on the ability to calibrate models to clinical data, efficiently and robustly. The current, most advanced 3D multi-physics models of cardiac functions incorporate representations of electrophysiology (EP), mechanics, and haemodynamics, which cover multiple scales—from the protein up to the organ scale [6,7,8,9,10,11,12]. However, in most current applications, these are, at least from a functional perspective, rather generic, and do not account for the known inter-individual variability among patients [13]. While such one-heart-fits-all modelling approaches are useful for gaining generic mechanistic insights, they are limited in their ability to diagnose, stratify, or plan therapies on a case-by-case basis, as evidence of correspondence between model and physiological reality in a given patient is lacking. These limitations must be addressed by creating personalised models that are calibrated to clinical data. The development of techniques to create such personalised models, often referred to as digital twins [5,14], models that replicate all available observations with high accuracy, constitute a significant challenge as parameter spaces of biophysically detailed models are high-dimensional, and available clinical data are afflicted with significant observational uncertainty and residual variability [13]. Beyond the high-dimensional parameter space, many model parameters cannot be observed clinically and, thus, must be estimated indirectly, based on quantities that are clinical observable; a procedure that can become computationally expensive.
In general, calibration of 3D multi-physics models is typically attempted by calibrating, in a first step, individual physics independently, and subsequent re-calibration, to account for secondary effects due to bi-directional coupling. For models of cardiac electromechanics (EM), typically EP components are calibrated first to match activation and repolarisation, before mechanical components are calibrated to data on pressure, volume, motion, or strain. Significant improvements in the personalisation of model components have been made recently in terms of robustness, automation and fidelity, including EP [15,16,17], afterload [18], and passive mechanics [19,20] models. In the present study, we focus on the fully automated personalisation of active mechanics models that characterise the active mechanical behaviour through corresponding constitutive equations. These can either be described in terms of stress (active stress approach) or in terms of strain (active strain approach) [21], with the first being more commonly used. The constitutive equations of the active stress approach involve the cellular active stress that is generated in cardiomyocytes through the mechanisms of excitation–contraction coupling [22]. There are numerous models available in the literature that describe the evolution of cellular active stress with varying degrees of mechanistic detail. Purely phenomenological low-dimensional models are often preferred in clinical modelling studies. These are able to account for salient phenomena, such as length dependence underlying the Frank–Starling mechanism and are easier to constrain with the limited data available in the clinic. They are simply driven by electrical activation time to trigger the onset of contraction, whereas biophysically detailed models require space-varying intracellular calcium concentration ( [ C a 2 + ] i ) traces as input. For a review on active mechanics models, the interested reader is referred to Niederer et al. [3].
Several fully automated personalisation approaches have been published for purely phenomenological models [18,23,24,25]. Here, the number of parameters that have to be calibrated is typically small, hence, the personalisation is computationally tractable. However, an increasing number of studies, e.g., pharmacological applications [26], look explicitly at the excitation–contraction coupling mechanisms for which purely phenomenological models are not suitable. For this reason, Longobardi et al. [27] introduced a fully automated personalisation approach that can also be used for biophysically detailed models. They applied Bayesian history matching (BHM) based on Gaussian process regression (GPR) models that emulate the expensive 3D organ-scale simulations at low computational cost. While this ensures that the actual personalisation process remains inexpensive, the cost involved in the generation of training data are substantial, which poses a limitation.
The aim our study was to develop an alternative fully automated and computationally efficient approach for the personalisation of biophysically detailed active mechanics models. To this end, a two-step multi-fidelity solution is suggested in Section 2. In brief, the active mechanical behaviour in a given 3D EM model is represented by a purely phenomenological model (low-fidelity model, LFM), which is personalised at the organ scale by calibration to clinical cavity pressure data. Then, median traces of nodal cellular active stress, [ C a 2 + ] i , and fibre stretch, are generated and utilised to personalise the desired biophysically detailed model (high-fidelity model, HFM) at the cellular scale using an 0D EM model. The personalisation approach was tested for a cohort of seven patient cases for which 3D models of human LV EM have been built previously [18], see Section 2.2. To be considered viable, the personalisation approach must demonstrate sufficient robustness against clinical data uncertainty, which we show for our approach in Section 3.3. Finally, we discuss the results in Section 4 and conclude the paper with a brief summary Section 5. All required equations to implement the methods in a software framework, as well as supplementary parameter values, are given in the Appendix A.

2. Materials and Methods

2.1. Clinical Data

Clinical data of all ( N = 7 ) aortic coarctation (CoA) patients from the CARDIOPROOF cohort (NCT02591940), see [18], were used. The institutional Research Ethics Committee approved the study following the ethical guidelines of the 1975 Declaration of Helsinki. Written informed consent was attained from the patients or their guardians.
The data of each patient include anatomical 3D-whole-heart (3DWH) magnetic resonance imaging (MRI) scans, an LV volume trace obtained from short-axis cine MRI scans, and LV pressure traces of several beats obtained from invasive catheterisation. A detailed description of the data acquisition process and clinical protocols used in this study are reported in [28]. To account for inconsistencies in the LV pressure and volume traces, a three-step pre-processing was performed previously [18]: firstly, the LV volume trace was adjusted to match LV volume data points that were derived from 3DWH-MRI scans acquired during diastasis; secondly, the LV pressure traces of all beats were averaged; thirdly, the averaged LV pressure trace was synchronised with the smoothed LV volume trace. The active mechanics personalisation approach is based on biomarkers and corresponding biomarker values of the LV pressure and volume trace were obtained as described in Appendix C.

2.2. 3d Model of Human Left Ventricular Electromechanics

2.2.1. Anatomical Model

A solid model of the LV and the aortic root (AR) was created based on 3DWH-MRI scans that were acquired during diastasis. Classification tags were applied to allow for local assignments of mechanical and electrical material properties and a finite element (FE) mesh was generated (Table 1). Since the model is not in an unloaded configuration, a backward displacement method [20,29] was applied first. Then, the unloaded configuration was inflated to the configuration at the clinical end-diastolic pressure and this was deemed the reference configuration [18]. To account for the unique fibre and sheet architecture of the LV, the principal eigenaxes along the fibre direction f 0 , the sheet direction s 0 , and the sheet-normal direction n 0 in the reference configuration were assigned to each element of the LV mesh using a previously developed Laplace–Dirichlet rule-based method [30]. Details of the semi-automated workflow for building anatomical models are described in [31,32].

2.2.2. Mechanical Model

The LV and the connected AR are represented by a deformable body B that consists of particles in the configuration Ω R 3 with Lipschitz boundary Ω . The deformation from the reference configuration Ω 0 ( X ) to the current configuration Ω t ( x ) is described by the deformation gradient F : = Grad x and the Jacobian J : = det F describes the change of the body’s volume. Furthermore, the right Cauchy-Green tensor is introduced to be C = F F . The tissue of the LV and the AR is assumed to be hyperelastic and nearly incompressible. To account for nearly incompressible behaviour, the deformation gradient is multiplicatively decomposed into the volumetric ( F vol ) and isochoric ( F ¯ ) parts [33]
F vol = J 1 / 3 I a n d F ¯ = J 1 / 3 F .
The spatiotemporal evolution of the displacement U in the LV and the AR is governed by Cauchy’s first equation of motion that reads
ρ 0 d 2 U d t 2 Div ( F S ) = 0 i n Ω 0 × ( 0 , T ) ,
with initial conditions
U = 0 a n d d U d t = 0 i n Ω 0 × { 0 } ,
where S is the second Piola–Kirchhoff stress, ρ 0 is the density in the reference configuration, d 2 U d t 2 is the acceleration, and d U d t is the velocity. To enforce physiological motion and to account for the cavity pressure that acts onto the endocardial surface, Robin spring boundary conditions and Neumann boundary conditions (Figure 1) are imposed on the boundaries Γ R 0 and Γ N 0 , respectively, with Γ R 0 ¯ Γ N 0 ¯ = Ω 0 ¯ .
More specifically, Robin boundary conditions refer to omni-directional springs that are applied to the aortic rim, and to uni-directional springs that are applied to the septum and to the epicardium. Latter mimics the effect of the pericardium that restricts changes of the outer shape of the heart. The pressure–volume relationship in the LV is described as in [34] and afterload is accounted for through an 0D-lumped three-element Windkessel model [35] (Figure 1) that describes the relationship between pressure and flow in the arterial system. This is coupled to the mechanical model during the ejection phase and patient-specific parameter values were already determined in previous work [18]. The relationship between pressure and flow across the aortic and mitral valve is described through a 0D diode representations (Figure 1) with forward resistances R AVf and R MVf , and backward resistances R AVb and R MVb , respectively. The backward resistances were set to 1000 mmHg ms mL 1 to prevent backflow.
Stress in the AR is assumed to arise only from the passive mechanical behaviour of the aortic tissue, i.e., S = S p . The passive mechanical behaviour is modelled by
S p = 2 Ψ C ,
where
Ψ = Ψ vol ( J ) + Ψ iso ( C ¯ )
is the strain energy function, where C ¯ is the isochoric right Cauchy–Green tensor. The strain energy function is additively decomposed into the volumetric contribution Ψ vol = κ 2 ( J 1 ) 2 with the penalty parameter κ and the isochoric contribution Ψ iso . AR tissue is assumed to be isotropic and is specified by the simple neo-Hookean strain energy function.
Stress in the LV arises not only from the passive but also from the active mechanical behaviour of ventricular tissue. Motivated by Hill’s muscle model passive and active stress are added (active stress approach) to obtain
S = S a + S p .
The passive mechanical behaviour is modelled by Equations (4) and (5) and the volumetric contribution is Ψ vol = κ 2 ln ( J ) 2 . LV tissue is assumed to be mechanically transversely isotropic and the isochoric contribution Ψ iso is specified by the strain energy function proposed by Guccione et al. [36]. Patient-specific values for the stiffness parameter C GUC were already determined in previous work [18]. Active stress is assumed to occur only in fibre direction f 0 and the active mechanical behaviour is thus modelled by
S a = S a ( f 0 · C f 0 ) 1 f 0 f 0 ,
where S a is the scalar cellular active stress. This is coupled to an 0D model describing the cellular active stress evolution. Here, the purely phenomenological Tanh model [37,38] and the biophysically detailed Land model [39] are used; detailed descriptions are given in Appendix A.1 and Appendix A.2.

2.2.3. Electrophysiological Model

The tissue of the LV is assumed to be electrically orthotropic and the spatiotemporal evolution of the transmembrane potential V m is governed by the reaction-eikonal (R-E) model [40]. The eikonal equation models the spatiotemporal evolution of electrical activation (wavefronts) and reads
( Grad t a ) v ( Grad t a ) = 1 in Ω 0 ,
with initial conditions
t a = t 0 in Γ 0 Ω 0 ,
where v is the (wavefront) velocity tensor and t a is a positive function that gives the electrical activation (wavefront arrival) times at the locations X Ω 0 . Electrical activation is initiated at the locations X Γ 0 in the vicinity of the septal, anterior, and posterior fascicles [34]. The development of the action potential upon electrical activation is modelled by the monodomain reaction–diffusion (R–D) equation    
β C m V m t = Div ( σ Grad V m ) + I foot ( t a ) β I ion in Ω 0 × ( 0 , T ) ,
with initial conditions
V m = V m res res in Ω 0 × { 0 } .
The electrical activation time determines the onset of a foot current density I foot , which mimics a subthreshold electrotonic current density that initiates the action potential starting out from the resting potential V m res . Velocities and conductivities σ along the three eigenaxes were set in accordance with [7].
The total ionic current density I ion is coupled to a 0D model that describes cellular EP. This usually also integrates a model of the [ C a 2 + ] i evolution that provides [ C a 2 + ] i as input for biophysically detailed models of cellular active stress evolution (Appendix A.2). The mammalian ventricular cardiomyocytes model of Luo and Rudy [41], from now on referred to as the LR1 model, was used. To produce a [ C a 2 + ] i trace that is in line with available experimental measurements in human cardiomyocytes, the purely phenomenological model of Rice et al. [42] was used and the parameter values were estimated based on experimental data reported in [43]. This model is from now on referred to as the Rice model and detailed descriptions including the parameter estimation strategy are given in Appendix A.4.
EP and mechanics are linked through various feedforward and feedback loops [44,45]. The action potential triggers active stress evolution through an increase of [ C a 2 + ] i (electromechanical coupling) and the resulting deformation modifies the transmembrane potential (mechano-electrical coupling). The latter can be caused through stretch-activated and stretch-modified ion currents and by stretch-induced changes of [ C a 2 + ] i . Stretch-induced changes of [ C a 2 + ] i modify not only the transmembrane potential but also the active stress (mechano-mechanical coupling). This arises from a modified affinity of troponin C for Ca2+ and from a modified sensitivity of the myofilaments to Ca2+. Electro- and mechano-mechanical coupling were considered but for the sake of simplicity, mechano-electrical coupling was not.

2.2.4. Spatiotemporal Discretisation and Numerical Solution

The spatial discretisation of the mechanical and EP model equations was done based on the Galerkin finite element (FE) method using tetrahedral elements. The same average mesh resolution of approximately 690 μ m was chosen for the two physics although EP processes are governed by smaller spatial scales than mechanical processes. This was made possible by the use of the R-E model that allows using a much coarser mesh resolution without causing numerical conduction slowing or block as observed in standard R-D models [40]. However, since EP processes are also governed by smaller temporal scales than mechanical processes, a different temporal discretisation was chosen. The time step for solving the EP and mechanical model equations were 0.025   m s and 1.0   m s , respectively. This choice of time step size is motivated by previous works using the same meshes [18]. To increase stability we considered Rayleigh damping of the generalised- α scheme used for the time integration of Equation (2), see [11]. Additionally, we made use of an approach developed by Regazzoni and Quarteroni [46] to stabilise velocity-dependent active stress models of cardiac mechanics. The cellular models always started out from their steady state computed for the clinical cardiac cycle length (duration of the clinical pressure trace; (Table 1)). If parameters were modified, the novel steady state was computed immediately. The spatial and temporal discretisation of the model and the solution of the arising systems of equations were realised in the FE framework Cardiac Arrhythmia Research Package (CARPentry) [40,47], built upon extensions of the openCARP EP framework [48] (http://www.opencarp.org, accessed on 30 January 2021). Numerical details for solving the EP [40,47,49] and the mechanical model equations [7] have been described in detail elsewhere. Both the solver components of the EP and the mechanical model have been verified in N-version benchmark studies [50,51]. Simulations were run on ARCHER2 (UK Research and Innovation) using 128 cores. The temporal discretisation of the cellular models and the solution of the arising systems of ordinary differential equations were realised with the tool bench included in openCARP. Simulations were run on a regular desktop computer using one core.

2.3. Active Mechanics Personalisation Approach

2.3.1. Step 1: Low-Fidelity Model Personalisation at the Organ Scale

In the first step, the active mechanical behaviour in the 3D EM models is described by a LFM. This is personalised at the organ scale by minimising the difference between simulated (sim) and clinical (clin) pressure biomarker values that is obtained from available clinical cavity pressure traces. Our specific choice of the LFM is the Tanh model [37,38]. The four major parameters maximum active stress, the time constants of the contraction (rise) and relaxation (decay) phase (rise time constant), and the duration of the cellular active stress transient are active stress biomarkers that are related to corresponding pressure biomarkers. Then, the personalisation can be performed by solving the following unconstrained minimisation problem:
min p LFM j = 1 4 B P clin , j B P sim , j ( p LFM ) B P clin , j 2 ,
where p LFM = { S max ref , τ S R ref , τ SD , t CR } are the parameters of the Tanh model (Appendix A.1) and B P = { P max , d P d t | max , d P d t | min , P T D 90 } are the pressure biomarkers (Appendix C) that are widely used, e.g., [18,26,27,37,38]). The formulation of the minimisation problems in this study is based on relative least squares as selected biomarkers vary significantly in magnitude; this would introduce an unintentional weighting of the terms otherwise. Owing to the relationship between parameters and biomarkers, a fixed point approach can be used to solve the problem with few iterations at low cost. The following choice of updates proved to be most promising for our application:    
S max ref i + 1 = S max ref i · P max clin P max sim , i , τ SR ref i + 1 = τ SR ref i · d P sim , i d t | max d P clin d t | max , τ SD i + 1 = τ SD i · d P sim , i d t | min d P clin d t | min , t CR i + 1 = t CR i · P T D 90 clin P T D 90 sim , i .
If, additionally, a clinical cavity volume trace is available, the personalisation can be extended to valve flow models (VFM) that describe the relationship between pressure and flow across valves. In this case, the minimisation problem reads
min p LFM , p VFM j = 1 4 B P clin , j B P sim , j ( p LFM ) B P clin , j 2 + k = 1 n VFM B V clin , k B V sim , k ( p VFM ) B V clin , k 2
where p VFM are the parameters of the VFM and B V are related volume biomarkers. Owing to the relationship between parameters and biomarkers, a fixed point approach can be used again for solving the problem. In the simple case of a diode valve with forward resistance R f , the additional update is:    
R f i + 1 = R f i · τ V clin τ V sim , i ,
where B V = τ V is the time constant of either the ejection (decay time constant τ VD ) or the filling phase (rise time constant τ VR ), see Appendix C.
A simulation with the current set of parameter values is performed in each iteration i until the minimisation problem is considered to be converged, i.e., the cost function reached a predefined threshold. Subsequently, traces of cellular active stress, [ C a 2 + ] i , and fibre stretch ( λ = f 0 · C f 0 ) are extracted for each node of the LV FE grid. These traces are first aligned by the respective nodal electrical activation times and then used to generate median traces S a ( t ) , [ C a 2 + ] i ( t ) , and λ ( t ) . Here, the median was chosen to be more robust against outlier traces.
In this study, the Tanh model and both the aortic and mitral VFM were personalised. These were integrated in the patient-specific 3D model of human LV EM (Section 2.2) and the required pressure and volume biomarker values were obtained from available patient-specific LV pressure and volume traces (Section 2.1). Initial guesses of the Tanh model are given in Table A1 and the forward resistances of the VFMs were initialised with 0.01 mmHg mL s 1 . One heart beat was simulated in each iteration and the threshold value for convergence of the minimisation problem was set to 0.1.

2.3.2. Step 2: High-Fidelity Model Personalisation at the Cell Scale

In the second step, the median nodal traces are used to personalise the desired 0D HFM at the cell scale. Biophysically detailed models of cellular active stress evolution are functions of [ C a 2 + ] i , the fibre stretch, and the fibre stretch rate. In theory, the personalisation could be performed by using the cellular active stress trace as target and the [ C a 2 + ] i trace, the fibre stretch trace, and the time derivative of the fibre stretch trace as input. However, the distribution of nodal fibre stretch in the myocardium at a given time point is very heterogeneous, which causes large errors in the personalisation. An exemplary distribution of nodal fibre stretch traces is shown in Figure 2. To overcome this issue, we suggest using an 0D EM models that is able to simulate the stretch of the cardiomyocyte, λ , during the cardiac cycle based on the equilibrium of active ( S a ) and passive stress ( S p ):
S a ( [ C a 2 + ] i , λ , d λ d t , t ) + S p ( λ , t ) = 0 f o r t ( 0 , T ) ,
with initial conditions
λ ( 0 ) = λ 0 .
The cellular active stress evolution is described by the HFM and the equilibrium equation can be solved in line with [46] to obtain the stretch trace and the stretch rate trace as time derivative. Here, the initial stretch λ 0 was set to the initial value of the median nodal trace and for an accurate personalisation of the active mechanics model it was found to be sufficient to match the minima λ min . To describe the cellular passive stress evolution, we used the three-element model of Land et al. [39] (Appendix A.3). Then, the minimum stretch can be controlled solely by the stiffness parameter a p . Please note that a p may be far from physiological when interpreted in the context of single cell stiffness. This is because fibre stretch traces produced at the organ scale are not only functions of local active and passive stress but also of the existing boundary conditions (Section 2.2.2) and deformations in the vicinity (Section 2.2.3). The only purpose of a p here is to produce a minimum stretch in line with the median nodal trace.
The personalisation can be performed by solving the following constrained minimisation problem:
min p HFM , a p j = 1 n HFM w S a j B S a tar , j B S a sim , j ( p HFM , a p ) B S a tar , j 2 + w λ B λ tar B λ sim ( p HFM , a p ) B λ tar 2 s . t . p min p p max , S a res = 0 , 0 S a ( t ) , d S a ( t ) d t | t > t S a max 0
where S a max = argmax t [ 0 , T ] S a ( t ) , p HFM are the parameters of the HFM, B S a are appropriate active stress biomarkers, and B λ = λ min is the stretch biomarker. The target biomarkers are obtained from the corresponding median nodal traces. In contrast to the first step, constraints are required since parameters p = { p HFM , a p } and the traces of cellular active stress and stretch can easily become non-physiological. The minimum requirements for the cellular active stress trace is that the resting stress is zero; that all stresses are equal or above zero; and that the stress rate is not positive after the maximum stress has been reached. If these requirements are met, no further constraints are required for the stretch trace. Non-negative weights w S a j and w λ were found to be helpful to control the impact of certain terms in order to obtain a physiological outcome.
Instead of using the median trace of nodal [ C a 2 + ] i as input for the HFM, the HFM can also be coupled to a model of [ C a 2 + ] i evolution that is integrated in the chosen cellular EP model. This comes with the advantage that the model of [ C a 2 + ] i  evolution can be included in the personalisation process. In this case, the constrained minimisation problem is extended to:
min p HFM , a p , p CEM j = 1 n HFM w S a j B S a tar , j B S a sim , j ( p HFM , a p , p CEM ) B S a tar , j 2 + w λ B λ tar B λ sim ( p HFM , a p , p CEM ) B λ tar 2 , s . t . p min p p max , S a res = 0 , 0 S a ( t ) , d S a ( t ) d t | t > t S a max 0 , 120   m s C T D 50 420   m s , 220   m s C T D 90 785   m s ,
where p CEM are the parameters of the [ C a 2 + ] i  evolution model. In this case, above mentioned constraints on the parameters p = { p HFM , p CEM , a p } and the cellular active stress trace were extended by constraints on the [ C a 2 + ] i trace. More specifically, C T D 50 and a C T D 90 —the [ C a 2 + ] i transient durations at the time of 50% and 90% decay from the maximum, respectively, measured from the electrical activation time—were required to be within the range given in [43,52]. 0D cell-scale simulations are computationally much less expensive than 3D organ-scale simulations which reduced the computational cost of the personalisation approach substantially. Further improvement of computational efficiency can be achieved by reducing the number of model parameters by means of a global sensitivity analysis (GSA), see Section 2.4.
In this study, the 0D model of cellular EM consisted of the Land models of cellular active stress (Appendix A.2) and passive stress evolution (Appendix A.3), and the Rice model that describes the evolution of  [ C a 2 + ] i on a purely phenomenological basis (Appendix A.4). Both the Land model of cellular active stress evolution and the Rice model were personalised and to this end, the following active stress biomarkers were used to set up the minimisation problem: B S a = { S a max , d S a d t | max , S T D 30 , S T D 50 , S T D 90 } (Appendix C).
In addition to the stiffness parameter a p , the nine most influential parameters that were identified by means of a GSA (see Section 2.4), were considered for estimation. Initial guesses and constraints on the parameters are given in Table A2, Table A3, Table A4 and Table A5 and the weights were w S a = 5 , 1 , 5 , 1 , 1 and w λ = 100 .
The constrained minimisation problem was solved as series of unconstrained problems by application of the penalty method (Appendix D). For solving, the population-based differential evolution method was used. This is a meta-heuristic global optimisation method [53] that was developed for multi-dimensional real-valued functions. Several studies [54,55] have demonstrated fast convergence, robustness, and good performance in real-world problems. As the gradient of the problem is not required, it can also be used for noncontinuous problems. This is advantageous for our application since some parameter value combinations within the admissible range may lead to failure of the 0D cell-scale simulation. The implementation of the differential evolution method was based on lmfit: non-linear least-squares minimisation and curve-fitting for Python [56] with default settings. In more detail, a Latin hypercube sampling was used to generate an initial population of candidate solutions that was then iteratively modified using the best1bin variant until convergence.
The two-step multi-fidelity approach is conceptually illustrated in Figure 2 and the general algorithm is given in Algorithm 1.

2.4. Global Sensitivity Analysis

The variance-based Sobol’ GSA [57,58] was used to identify the parameters that the chosen active stress biomarkers are most sensitive to in the 0D EM model. In general, when the model inputs are varied, it measures the sensitivity of the model outputs to each model input by the fraction of variance attributed to each model input alone (first-order sensitivity index S 1 ) or by the fraction of variance attributed to each model input including interactions with the other model inputs (total-effect sensitivity index S T ). The Sobol’ GSA is based on an all-at-a-time sampling strategy and the sensitivity indices range from 0 (no sensitivity) to 1 (maximum sensitivity).
Algorithm 1 Two-step multi-fidelity approach for personalising biophysically detailed active mechanics models.
1:
Initialise counter i = 0 and parameters of the low-fidelity model p LFM 0 , see Appendix A.1.
2:
do
3:
    Solve 3D organ-scale model with low-fidelity model and parameters p LFM i .
4:
    Compute error between simulated ( B p sim ) and clinically measured ( B p clin ) pressure biomarkers.
5:
    Update parameter set using fixed point approach, see Equations (13) and (15), to get p LFM i + 1 .
6:
    Update counter i = i + 1 .
7:
while error > threshold  
8:
Generate median nodal traces of active stress ( S a ( t ) ), intracellular calcium concentration
( [ C a 2 + ] i ( t ) ), and fibre stretch ( λ ( t ) ) from low-fidelity model results.
9:
Compute target biomarkers for active stress ( B S a tar ) and fibre stretch ( B λ tar ) from S a ( t ) and λ ( t ) ,
respectively.  
10:
Initialise counter j = 0 and parameters of the high-fidelity model p HFM 0 , see Appendix A.2.
11:
do
12:
    Solve low-cost 0D cell-scale model with high-fidelity model using p HFM j and [ C a 2 + ] i ( t ) .
13:
    Compute error between simulated ( B S a sim , B λ sim ) and target ( B S a tar , B λ tar ) biomarkers.
14:
    Update parameter set based on differential evolution method to get p HFM j + 1 .
15:
    Update counter j = j + 1 .
16:
while error > threshold
In this study, the Sobol’ GSA was performed to identify the nine parameters of the Land and the Rice model (inputs) that the five active stress biomarkers (Section 2.3.2; outputs) are most sensitive to in the 0D EM model. Saltelli’s sampling scheme [59] with N = 1024 was applied to generate 45,506 parameter samples. Lower and upper bounds of the parameters were in line with the parameter constraints in the second step of the personalisation approach (Section 2.3.2). The simulations were performed with a cardiac cycle length of 1000 ms and to produce stretch traces in line with those seen in the organ-scale simulations, the initial stretch λ 0 was set to 1.1 and the stiffness factor a p was set to 42 kPa leading to comparable minima.
Samples that produced non-physiological cellular active stress traces were excluded. The exclusion criteria were adopted from the cellular active stress constraints in the second step of the personalisation approach (Section 2.3.2) with the only difference that the resting stress was allowed to take values up to 10% of S a max to increase the number of data points for the analysis. Since the number of output and input values must be equal, the biomarkers of the non-physiological traces were set to the means of the biomarkers of the physiological traces.
Finally, the parameters were ranked by their sensitivity. To this end, the total-effect sensitivity indices with respect to each biomarker were added. The implementation of the GSA was based on SALib—Sensitivity Analysis Library in Python [60].

2.5. Robustness Analyses

2.5.1. Clinical Data Uncertainty

To quantify the robustness of the active mechanics personalisation approach against uncertainties that clinical data are afflicted with, their effect on the estimated parameter values and consequences for the simulated pressure and volume traces were analysed (uncertainty propagation). First, samples of clinical biomarker values were generated that evenly filled a range of ±10% around the measured values. Then, the samples were used to perform the active mechanics personalisation.

2.5.2. Initial Guess Variation

To quantify the robustness of the active mechanics personalisation approach against variation of initial guesses, their effects on the estimated parameter values, and consequences for simulated pressure and volume traces were analysed. For this purpose, samples of initial guesses that evenly filled a range of ±50% around the default values were generated. This was done for both steps.
The sampling was performed using Latin hypercube sampling implemented in SALib [60]. Owing to the marked computational costs involved, the robustness analyses were carried out only for the case 02-CoA and a sample size of ten and five, respectively.

3. Results

3.1. Global Sensitivity Analysis

Figure 3 shows the sensitivity of the chosen active stress biomarkers to the parameters of the Land and the Rice model in the 0D EM model. The comparison of first-order ( S 1 ) and total-effect indices ( S T ) demonstrate that the influence of the parameters on the active stress biomarkers was primarily caused by interactions among them. For the given admissible parameter ranges, the Land model parameters [ C a 2 + ] 50 ref and n TRPN were most influential across all biomarkers, followed by the Rice model parameters [ C a 2 + ] res and [ C a 2 + ] max . This highlights the role of troponin C kinetics and the [ C a 2 + ] i dynamics in the evolution of cellular active stress. The role of the [ C a 2 + ] i dynamics is further emphasised by the fact that all parameter of the Rice model were among the nine most influential parameters: τ CD was ranked sixth and τ CR was ranked eight. The list was completed by the Land parameters β 0 , n Tm , and k UW . These nine parameters were considered in the second step of the active personalisation approach.

3.2. Patient Cohort Results

The active mechanics personalisation approach, as described in Section 2.3, was applied to all N = 7 patient cases: 01-CoA–07-CoA. In the first step, LV pressure traces were used to personalise the Tanh model and LV volume traces were used to include the aortic and the mitral VFM in the personalisation process. This was done based on the 3D model of human LV EM. For six out of seven cases, two to four iterations were needed for convergence (Table 1). For case 04-CoA, the simulation aborted at the tenth iteration and therefore, the results of the ninth iteration were taken. The convergence behaviour for this particular case is shown in Figure A2. The estimated parameter values for all models are given in Table A6.
Figure 4 compares simulated and clinical LV pressure–volume loops and the individual LV pressure and volume traces. Moreover, Table 2 compares simulated and clinical values of relevant pressure and volume biomarkers ( P max , d P d t | max , d P d t | min , P T D 90 , S V ). The personalised models were able to reproduce the clinical pressure with good agreement. The relative differences between the simulated and clinical values of P max and P T D 90 had means (standard deviations; SD) of 5.4% (SD: 3.4%) and 1.4% (SD: 2.2%), respectively. For d P d t | max and d P d t | min , the mean relative differences were slightly larger with 13.6% (SD: 8.0%) and 7.7% (SD: 8.0%), respectively. The brief pressure drop after the first peak seen in the clinical data are considered a measurement artefact. Except for case 04-CoA, the personalised models were also able to reproduce the clinical volume with good agreement, albeit the agreement in the systolic phase was better than in the diastolic phase. The mean relative difference between the simulated and clinical S V was 10.0% (SD: 5.9%). However, simulated end-diastolic volumes were substantially smaller than clinical end-diastolic volumes. Except for case 04-CoA, the good agreement between simulated and clinical pressure and volume translated into a good agreement in the pressure–volume relationship.
Median traces of nodal cellular active stress and fibre stretch obtained from the simulations of the personalised 3D models of human LV EM were then used as targets to personalise the Land model and the Rice model in the second step. The personalisation was done based on the 0D model of cellular EM and together with the stiffness parameter, a total number of ten parameters were included in the fitting. The number of iterations until convergence was between 8561 and 20,766 (Table 1). The estimated parameters are given in Table A7. Figure 5a compares the median nodal cellular active stress traces obtained in the first step with the cellular active stress traces obtained in the second step. For all cases, the maximum active stresses and transient durations at 30% decay from the maximum were in line with each other. This is a consequence of the weighting in the cost function. In contrast, the shapes in the plateau phase were slightly different and the resting value was achieved later after personalisation in the second step. Figure 5b compares the [ [ C a 2 + ] i traces. Please note that the median traces of nodal [ C a 2 + ] i are not targets for the personalisation in the second step. For all cases, the resting value after personalisation in the second step was at the lower bound and the maximum value was at the upper bound. This allowed larger transient durations of the cellular active stress traces. The stretch traces are not shown since the only important information is whether the minimum could be matched. Minimum fibre stretches of the median nodal traces ranged from 0.82 to 0.86 and matches were found for all cases. To this end, the stiffness parameter a p of the Land model of cellular passive stress evolution was fitted between 47.3 kPa and 86.5 kPa with a mean of 62.4 kPa (SD: 12.5 kPa).
Finally, the personalised models were incorporated in the 3D EM model to simulate one heart beat and to compare pressure and volume traces. The simulated pressure–volume loops and the individual pressure and volume traces were in good agreement with those produced by the personalised models in the first step and consequently also with the clinical data (Figure 4). The mean relative differences between the simulated and clinical values of P max and P T D 90 remained almost unchanged: 5.5% (SD: 2.3%) and 2.9% (SD: 1.5%), respectively. For d P d t | max and d P d t | min , the mean relative differences increased 22.5% (SD: 21.8%) and 39.9% (SD: 21.3%), respectively. The mean relative difference between the simulated and clinical S V also remained almost unchanged 11.1% (SD: 6.2%).

3.3. Robustness Analyses

3.3.1. Clinical Data Uncertainty

The robustness against clinical data uncertainty was analysed for the case 02-CoA. For this purpose, the active mechanics personalisation was performed based on ten samples of clinical biomarker values that were within a range of ±10% around the measured values. Table 3, Table 4, Table 5 and Table 6 give the estimated parameter values, the resulting pressure and volume biomarker values, and the respective differences relative to the original values for the first and the second step. In the first step, the mean relative differences of the Tanh model parameter values were all below 8.0% and the mean relative differences of the VFM parameter values were similar with 10.8% (SD: 8.1%) for R AVf and 5.0% (SD: 3.3%) for R MVf . The relative differences of the resulting pressure and volume biomarker values were all below 6.5%. In the second step, the mean relative difference for the Land parameters [Ca 2 + ] 50 ref and β 1 was comparable to the mean relative difference for the Tanh parameters: 1.3% (SD: 3.3%) and 0.3% (SD: 0.2%). For the Land model parameters k UW , n TRPN , n Tm , the mean relative difference was larger: 32.7% (SD: 40.4%), 27.0% (SD: 16.4%), 11.5% (SD: 12.4%). The mean relative differences of the Rice parameter values and the stiffness parameter a p (results not shown) were all below 5.4%. The relative differences of the resulting pressure and volume biomarker values were below 7.0% except for d P d t | max (10.5%; SD: 7.3%). The pressure–volume loops and the individual pressure and volume traces are compared in Figure 6.

3.3.2. Initial Guess Variation

The robustness against variation of initial guesses was analysed for case 02-CoA. Five samples of initial guesses within a range of ±50% around the default values were used to perform the active mechanics personalisation. Table 7, Table 8, Table 9 and Table 10 give the estimated parameter values, the resulting pressure and volume biomarker values, and the respective differences relative to the results produced by the default initial guesses for the first and the second step.
In the first step, the mean relative differences for the Tanh model parameters S m a x ref , τ SD , T CR , and the VFM parameter R MVf were below 2.9%. For τ S R ref (10.9%; SD: 11.5%) and in particular for R AVf (30.4%; SD: 28.5%), they were larger.
In the second step, the mean relative differences for the Land model parameters β 1 and [Ca 2 + ] 50 ref , for the Rice model parameters [ C a 2 + ] res , [ C a 2 + ] max , and for the stiffness parameter a p (results not shown) were below 2.8%. For the remaining half of the parameter set, they were up to 59.3%. Large relative differences were mainly due to outliers in individual samples which is also indicated by the corresponding large standard deviations (up to 85.2%). The mean relative differences of the resulting pressure and volume biomarker values were below 3.4% in the first step and below 12.8% in the second step.

4. Discussion

This study describes a novel approach for the fully automated personalisation of biophysically detailed active mechanics models in 3D EM models. Motivated by the aim to keep the computational cost low, we suggest a two-step multi-fidelity solution based on clinical cavity pressure data. Here, the purely phenomenological Tanh model [37,38] was used as LFM and the biophysically detailed Land model [39] was used as HFM. The personalisation approach was tested for seven patient cases with previously built 3D models of human LV EM. These account for all cases in the CARDIOPROOF cohort presented by Marx et al. [18] for which invasively measured LV pressure data were recorded. Since volume traces were also available, the personalisation of the Tanh model in the first step was extended to an aortic and a mitral VFM. The personalisation of the Land model in the second step was extended to the Rice model [42] that represents [ C a 2 + ] i evolution on a purely phenomenological basis. This was done because information on [ C a 2 + ] i traces in vivo are missing but can have considerable influence on the parameter estimation in biophysically detailed models as demonstrated in Tøndel et al. [61].

4.1. Computational Cost

For the six successfully converged patient cases, only two to four 3D organ-scale simulations were required in the first step to personalise the Tanh model and the aortic and a mitral VFM (Table 1). The Land and the Rice model were personalised based on computationally much less expensive, single-core 0D cell-scale simulations. Here, we achieved a reduction of computational costs by reducing the number of considered model parameters to those that the cellular active stress trace is most sensitive to. To this end, a GSA was performed which demonstrated that the active stress is most sensitive to parameters related to the troponin C kinetics and [ C a 2 + ] i dynamics (Figure 3). This is in line with previous observations by Tøndel et al. [61]. Overall, ten parameters were included in the fitting and the number of 0D cell-scale simulations was between 8561 and 20,766 (Table 1).

4.2. Goodness of Fit and Robustness

For the six successfully converged patient cases, the agreement between simulated and clinical LV pressure and volume traces was good for both the personalised Tanh model and the personalised Land model (Figure 4, Table 2). Substantial differences were only found for end-diastolic volumes; this was to be expected because the 3D model of human LV EM does not account for the atrial kick. This limitation could be addressed by coupling the 3D model of human LV or biventricular (biV) EM to a closed-loop lumped 0D model, see, e.g., [11,12], which represents the function of the remaining chambers and the circulation. However, this would require estimation of additional parameters of the closed-loop model which is beyond the scope of this work.
Mirams et al. [13] highlighted the importance of considering uncertainty in real-world data when calibrating models. In specific, observational uncertainty accounts for errors in the measurement process and residual uncertainty describes intrinsic and extrinsic variability in the biological system that is studied. We therefore tested the robustness of the presented personalisation approach against clinical data uncertainty which was estimated to be around ±10% of the measured values in line with previous work [18,20]. While variations in the estimated parameter values and resulting biomarkers are to be expected, these should ideally be not larger than the input data variation. Indeed, the mean differences relative to the original values were below 8.0% for the Tanh model parameters, below 10.8% for the VFM parameters (Table 3), and below 5.4% for the Rice model parameters and the stiffness parameter a p . Mean relative differences of all but two Land model parameters (Table 5) were below 11.5%. The maximum was 32.7%. Moreover, the consequences for the pressure and volume biomarkers were small: mean differences relative to the original values were below 6.5% in the first step (Table 4) and below 10.5% in the second step of the personalisation approach (Table 6, Figure 6). Overall, this demonstrates a high robustness against clinical data uncertainty.
In addition, the robustness against variation of initial guesses (by ±50%) was tested. It was found to be high in the first step. Except for the Tanh model parameter τ S R ref (10.9%) and the VFM parameter R AVf (30.4%), the mean differences relative to the results produced by the default initial guesses were below 2.9% (Table 7) and the consequences for the pressure and volume biomarkers were small with mean differences below 3.8% (Table 8). The robustness was found to be lower in the second step. For two Land model parameters ( β 1 , [ C a 2 + ] 50 ref ), two Rice model parameters ( [ C a 2 + ] res , [ C a 2 + ] max ), and the stiffness parameter a p , the mean relative differences were below 2.8%, but they were up to 59.3% for the other parameters. However, large mean relative differences resulted mainly from outliers (Table 9). The consequences for the pressure and volume biomarkers were larger with mean differences up to 12.8% (Table 10). The robustness to initial value variation of the second step may be improved by considering other variants of the differential evolution method [53].

4.3. Comparison to Other Active Mechanics Personalisation Approaches

Kayvanpour et al. [23] presented an personalisation approach for an active mechanics model that was integrated in 3D models of human biV EM. They used a purely phenomenological model published in Sermesant et al. [62] to describe the evolution of cellular active stress. The parameter that controls the maximum was estimated together with a parameter of the passive stress based on clinical ejection fraction, stroke volume, end-diastolic and end-systolic volume, and end-diastolic and end-systolic pressure. Asner et al. [24] presented a personalisation approach for an active mechanics model that was integrated in 3D models of human LV mechanics and the cellular active stress was multiplicatively decomposed into a reference stress and a factor that accounts for stretch dependency. They treated the reference stress as parameter, which they estimated for various time points in the cardiac cycle based on clinical wall displacements and cavity pressures. In contrast to the wall displacements, the cavity pressures were not measured but derived from an empiric reference trace [63] that was adjusted based on estimates of the end-diastolic and the maximum pressure, and the time points of mitral and aortic valve opening and closing. Finsberg et al. [25] presented a personalisation approach for an active mechanics model that was integrated in 3D models of human biV mechanics. They estimated the local cellular active stress directly over time based on clinical cavity volumes and circumferential strains. Moreover, they presented an analogue for the active strain approach and estimated the local cellular active strain over time instead. Marx et al. [18] presented a personalisation approach similar to the first step of our approach. Integrated in the 3D models of human LV EM that were also used in this study, they estimated the parameters of the Tanh model [37,38] based on the corresponding clinical pressure biomarkers. However, in contrast to our approach, they set the duration of the cellular active stress transient to be the RT interval of an available clinical electrocardiogram because they were only interested in the systolic phase.
The vital difference between our approach and the described previous works is that it cannot only be used for purely phenomenological but also for biophysically detailed models. To the best of the authors knowledge, only one other approach for that purpose has been published. Longobardi et al. [27] applied BHM based on GPR models to personalise an active mechanics model that was integrated in a 3D model of rat biV EM. The evolution of cellular active stress was described by the biophysically detailed model published in Land et al. [64]. While GPR models are even cheaper to solve than 0D models, the generation of a sufficiently large training data set for highly accurate GPR models requires a high number of 3D organ-scale simulations. This limits the efficiency as the number of 3D organ-scale simulations predominantly determines the computational cost of the personalisation approach. (Considering eight parameters), Longobardi et al. [27] reported on 1024 simulations to generate training data for the GPR models being used in the first wave of BHM and at each successive wave, the GPR models were updated by training data that is generated by another 256 simulations. In comparison, only between two and four 3D organ-scale simulations were required in the presented approach.

4.4. Limitations

While the results of this study are very promising to further the development of cardiac digital twins, some limitations have to be considered. First, the clinical cavity pressure data used for the personalisation approach can only be collected invasively which not only entails considerable efforts but is also associated with certain risks for the patient. Nevertheless, as was done previously for the LV [24], the required data could also be derived by utilising an empiric reference pressure trace that is adjusted based on non-invasive measurements [63]. For patient-specific adjustments, [63] used the time points of aortic and mitral valve opening and closing determined by echocardiography and [24] further used estimates of the end-diastolic [65,66], and the maximum pressure [67] based on phase-contrast MRI and cuff pressure measurements, respectively. Second, the parameter constraints were selected based on two studies [39,61] or were set to ±50% of the original values. Constraining parameters in personalisation approaches is crucial and, therefore, it is desirable to extend the database in the future. Third, the approach was tested for a rather small cohort of seven patient cases and only for the active mechanics of the LV. Kayvanpour et al. [23], Finsberg et al. [25], Longobardi et al. [27] applied their personalisation approach to the active mechanics of both ventricles and an extension of our approach to multiple chambers could also be done straightforwardly. Yet, this work focuses on the methodology of the novel multi-fidelity approach and tests for larger cohorts and more chambers will be left to future studies. In this regard, machine learning techniques could be employed to handle the challenges that arise in big data [68] and, hence, to further speed up parameter calibration. Fourth, as data and results were used from a previous work by Marx et al. [18], we also applied active stresses only in fibre direction. However, several studies [69,70,71] have shown that active stresses in the cross-fibre direction can be as large as 40% of those in fibre directions due to a dispersion of cardiomyocytes in tissue. This was already considered in Finsberg et al. [25] and recently in a mechanistically more accurate mechanical framework in Augustin et al. [11] and we see no obstacle in using other active mechanical models for the presented approach. Fifth, the representation of cellular EP and the [ C a 2 + ] i evolution was simplified. There are numerous biophysically detailed models of cellular EP and [ C a 2 + ] i  evolution available, e.g., [72], and if mechanistic detail is crucial for the study purpose, the simplified models have to be substituted. Finally, although we could show that the approach is widely robust to variations of initial guesses, we cannot prove uniqueness of the estimated parameters. Further, we cannot provide a rigorous proof that the fixed point approach given in Section 2.3.1 will always converge. In one case, we also encountered parameter values that caused the 3D organ-scale simulation to fail. Here, parameter constraints may improve the stability.

5. Conclusions

A novel approach for the fully automated personalisation of biophysically detailed active mechanics models was presented. The great strength of this approach lies in its efficiency with only a few 3D organ-scale simulations needed. Further, the application to a cohort of seven patient cases demonstrated an accurate and robust estimation of model parameters that resulted in good agreement of simulated and clinical LV pressure and volume traces. This also holds true when uncertainty in the clinical data and variations of initial guesses were taken into account. Thus, the presented workflow constitutes a further step forward towards the personalised modelling of active mechanical behaviour in the heart. As such, this approach is considered highly suitable for integration in workflows for building digital twins of cardiac EM—from a single patient to an entire cohort.

Author Contributions

Conceptualisation, A.J., C.M.A., G.P.; methodology, A.J.; software, A.J., M.A.F.G., C.M.A.; validation, A.J.; formal analysis, A.J.; investigation, A.J.; data curation, A.J., M.A.F.G.; writing—original draft preparation, A.J.; writing—review and editing, A.J., M.A.F.G., C.M.A., G.P.; visualisation, A.J., M.A.F.G.; funding acquisition, A.J., C.M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received funding from the European Union’s Horizon 2020 research and innovation programme under the ERA-NET co-fund action no. 680969 (ERA-CVD SICVALVES) funded by the Austrian Science Fund (FWF), grant I 4652-B to C.M.A., from the Austrian Science Fund (FWF) grant I2760-B30, and from BioTechMed-Graz, Austria, as part of the BioTechMed-Graz Flagship Award ILearnHeart to G.P. Furthermore, A.J. was supported by the German Research Foundation (Walter Benjamin Fellowship, no. 468256475). Simulations for this study were performed on ARCHER2 (http://www.archer2.ac.uk/, accessed on 30 January 2022), the UK national high-performance computing service located at the University of Edinburgh using core hours received from the eCSE funding program (application ID ARCHER2-eCSE01-21). Open Access Funding by the Austrian Science Fund (FWF).

Institutional Review Board Statement

The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Ethics Committee of Charité—Universitätsmedizin Berlin (no. EA2/172/13).

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Acknowledgments

The authors thank Thomas Grandits for the helpful discussions on the constrained minimisation problem in the second step of the personalisation approach.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
3DWHthree-dimensional-whole-heart
ARaortic root
BHMBayesian history matching
biVbiventricular
CoAaortic coarctation
CEMintracellular calcium concentration evolution model
EMelectromechanics
EPelectrophysiology
FEfinite element
GSAglobal sensitivity analysis
GPRGaussian process regression
HFMhigh-fidelity model
LFMlow-fidelity model
LVleft ventricle/ventricular
MRImagnet resonance imaging
SDstandard deviation
VFMvalve flow model

Appendix A. Models of Cellular Mechanics and the Intracellular Calcium Concentration Evolution

Appendix A.1. Tanh Model of Cellular Active Stress Evolution

The Tanh model [37,38] describes the evolution of cellular active stress as function of fibre stretch and the electrical activation time:
S a ( t a , λ ) = 0 if t t s ( t a , t emd ) S m a x ref ϕ tanh 2 t τ SR tanh 2 t CR t τ SD if t s ( t a , t emd ) < t t e 0 if t t e .
Here, S m a x ref is the maximum isometric cellular active stress and
ϕ = tan a 6 ( λ a 7 )
is a nonlinear function that describes stretch effects on the generated active stress. The coefficient a 6 corresponds to the degree of the stretch dependency and the coefficient a 7 is the fibre stretch below which no active stress can be generated. Furthermore,
τ SR = τ SR ref + a 4 ( 1 ϕ ) .
is the time constant of the contraction phase (rise time constant) that accounts for stretch effects. The isometric value is τ SR ref and the coefficient a 4 corresponds to the degree of the stretch dependency. The time constant of the relaxation phase (decay time constant) is denoted by τ SR and t CR is the duration of the entire contraction–relaxation cycle (transient). The contraction–relaxation cycle starts at t s and ends at t e . The starting time is the electrical activation time t a plus some electromechanical delay t emd to account for the time lag between electrical activation and the onset of contraction. The parameter values are given in Table A1.
Table A1. Parameters of the Tanh model [37,38].
Table A1. Parameters of the Tanh model [37,38].
ParameterUnitValue
S a ref (kPa)100
τ SR (ms)40
τ SD (ms)110
t CR (ms)550
a 4 (ms)500
a 6 (-)5
a 7 (-)0.7
t emd (ms)15

Appendix A.2. Land Model of Cellular Active Stress Evolution

The Land model [39] describes the evolutionn of cellular active stress as function of [ C a 2 + ] i and both the fibre stretch λ and the fibre stretch rate d λ d t . It is composed of a model of thin filament kinetics and a model of the cross bridge cycle. The model of thin filament kinetics describes the interactions of Ca2+, troponin C, troponin I, and tropomyosin that control the availability of myosin binding sites on actin. The dynamics of the interaction between Ca2+ and troponin C is described by
d C a T R P N d t = k TRPN [ C a 2 + ] i [ C a 2 + ] 50 ( λ ) n TRPN ( 1 C a T R P N ) C a T R P N ,
where C a T R P N is the fraction of regulatory troponin C sites with bound Ca2+, k TRPN represents the unbinding rate of Ca2+ from troponin C, and n TRPN is the cooperativity of the binding between Ca2+ and troponin C. The value of [ C a 2 + ] i at which half of the maximum active stress generated is denoted by [ C a 2 + ] 50 . Stretch effects are phenomenologically captured by
[ C a 2 + ] 50 = [ C a 2 + ] 50 ref + β 1 ( min ( λ , 1.2 ) 1 ) ,
where [ C a 2 + ] 50 ref is the isometric value that is scaled by β 1 . The fraction of regulatory troponin C sites with bound Ca2+ drives the unblocking of tropomyosin:
d B d t = k B C a T R P N n Tm 2 U k U C a T R P N n Tm 2 B ,
where B is the fraction of blocked myosin binding sites on actin,
k B = k U C a T R P N n Tm 1 r S ( 1 r S ) r W
and k U are the troponin I and tropomyosin rate constants, respectively, n Tm is the steady-state relation between C a T R P N and the fraction of unblocked binding sites ( 1 B ), and U is the fraction of the unblocked myosin binding sites with no cross bridges formed.
The model of the cross bridge cycle accounts for three states: the unbound, the weak (pre-powerstroke), and the strong (post-powerstroke) state. It reads
U = ( 1 B ) S W , d W d t = k UW U k WU W k WS W γ WU W , d S d t = k WS W k SU S γ SU S ,
where W and S are the weak, and the strong states, respectively, and k UW , k WU , k WS , k SU and are transition rates. Latter are defined by
k WS = k UW 1 r W 1 k WS , k SU = k WS r W 1 r S 1 ,
with the steady-state ratios
r W = steady - state W U + W , r S = steady - state S U + W + S .
The distortion-depending unbinding rates of the cross bridges are
γ WU = γ W | ξ W | , γ SU = γ S ( ξ S 1 ) if ξ S + 1 < 0 γ S ξ S if ξ S + 1 > 1 0 otherwise ,
and these are coupled to a distortion-decay model given by
d ξ W d t = A W d λ d t c W ξ W , d ξ S d t = A S d λ d t c S ξ S .
Here, ξ W and ξ S are the stretch rate-dependent mean distortions and
A W = A S = A eff r S ( 1 r S ) r W + r S
are related to the magnitude of the instantaneous response to the distortion with some scaling A eff , whereas
c W = Φ k UW U W , c S = Φ k WS W S
are related to the magnitude of the decay rate of the distortion. The introduction of Φ eliminates a parameter by reducing two parameters ( c W , c S ) to one.
Finally, the active stress is given by
S a ( [ C a 2 + ] i , λ , d λ d t ) = h ( λ ) S m a x ref r s [ ( ξ S ( d λ d t ) + 1 ) S + ξ W ( d λ d t ) W ] ,
where S m a x ref is the maximum isometric active stress and fibre stretch effects on the generated active stress are phenomenologically captured by
h = 1 + β 0 ( λ + min ( λ , 0.87 ) 1.87 ) ,
where β 0 represents the change in maximum active stress based on changes in filament overlap. The parameter values are given in Table A2.
Table A2. Parameters of the Land model [39] (active).
Table A2. Parameters of the Land model [39] (active).
ParameterUnitValue
k TRPN (ms 1 )0.1
n TRPN (-)2
[ C a 2 + ] 50 ref ( μ M)0.805
k U (ms 1 )1
n Tm (-)5
T R P N 50 ()0.35
k UW (ms 1 )0.182
k WS (ms 1 )0.012
r W (-)0.5
r S (-)0.25
γ S (-)0.0085
γ W (-)0.615
Φ (-)2.23
A eff (-)25
β 0 (-)2.3
β 1 (-)−2.4
S m a x ref (kPa)120

Appendix A.3. Land Model of Cellular Passive Stress Evolution

The Land model [39] describes the evolution of cellular passive stress by a three-element model similar to a standard linear solid. It consists of an elastic spring (E1) in parallel to another elastic spring (E2) in series with a viscous dashpot (V) and reads
S p ( λ ) = S E 1 + S E 2 = S E 1 + S V ,
with the stress components
S E 1 = a p ( e p b C 1 ) , S E 2 = a p k p C S ,
S V = a p η l d C V d t if d C V d t > 0 a p η S d C V d t if d C V d t < 0 ,
the series strain constraint
C = C S + C V ,
and the series stress constraint
S E 2 = S V .
Strain C is defined as C = λ 1 with stretch defined as the ratio of the current and initial cardiomyocyte length. The stiffness parameter a p is included in all stress components, such that it is suitable for scaling the passive stress. The parameter values are given in Table A3.
Table A3. Parameters of the Land model [39] (passive).
Table A3. Parameters of the Land model [39] (passive).
ParameterUnitValue
a p (kPa)2.1
b p (-)9.1
k p (-)7
η l (ms 1 )200
η S (ms 1 )20

Appendix A.4. Rice Model of the Intracellular Calcium Concentration Evolution

The Rice model [42] describes the evolution of  [ C a 2 + ] i as a function of the electrical activation time by
[ C a 2 + ] i ( t a ) = [ C a 2 + ] res if t t a [ C a 2 + ] max [ C a 2 + ] res β exp t t s τ C R exp t t s τ C D + [ C a 2 + ] res if t > t a ,
with
β = τ CR τ CD 1 τ CR τ CD 1 τ CR τ CD 1 1 τ CR τ CD ,
where [ C a 2 + ] res and [ C a 2 + ] max are the resting and the maximum [ C a 2 + ] i , respectively, and τ CR and τ CR are the time constants during the rise and the decay phase of the [ C a 2 + ] i transient, respectively. Starting time of the transient is the electrical activation time: t s = t a . The model was calibrated based on an experimentally measured [ C a 2 + ] i trace in human cardiomyocytes [43]. To this end, the unconstrained minimisation problem
min p Rice j = 1 n [ C a 2 + ] e x p , j [ C a 2 + ] sim , j ( p Rice ) 2 ,
was solved to estimate the parameters p Rice . Here, j = 1 , , n are all data points of the trace. Powell’s method implemented in the library lmfit: Non-linear least-squares minimisation and curve-fitting Python library [56] was used. The estimated parameter values are given in Table A4.
Table A4. Parameters of the Rice model [42].
Table A4. Parameters of the Rice model [42].
ParameterUnitValue
[ C a 2 + ] res ( μ M)0.15
[ C a 2 + ] max ( μ M)0.6
τ CR (ms)129
τ CD (ms)128

Appendix B. Parameter Bounds

Physiological lower and upper bounds were applied to the parameters of the Land model of cellular active and passive stress evolution and the Rice model. The bounds for the parameters of the Land model of cellular active stress evolution were either set according to Tøndel et al. [61] or according to Land et al. [39] or they were set to ±50% of the original value if no information were available. The bounds for the Rice model parameters were set to ±50% of the original value in agreement with the [ C a 2 + ] i traces of the experimentally calibrated model cohort shown in Passini et al. [52]. The range of the stiffness factor in the Land model of cellular passive stress evolution was set to be very wide as this parameter has no physiological meaning in the context of the second step of the active mechanics personalisation approach. All bounds are given in Table A5.
Table A5. Parameter bounds used in this study. Lower bounds are denoted by LB and upper bounds are denoted by UB. If applicable, the source is given.
Table A5. Parameter bounds used in this study. Lower bounds are denoted by LB and upper bounds are denoted by UB. If applicable, the source is given.
ParameterUnitLBUBSource
Land model (active)
k TRPN (ms 1 )0.050.4[61]
n TRPN (-)15[61]
[ C a 2 + ] 50 ref ( μ M)0.52.0[61]
k u (ms 1 )0.012.00[39]
n Tm (-)37[39]
T R P N 50 (-)0.30.5[61]
k UW (ms 1 )0.0260.312[39]
k WS (ms 1 )0.0040.048[39]
r W (-)0.250.75[39]
r S (-)0.100.25[39]
γ S (-)0.0050.020[39]
γ W (-)0.053.00[39]
Φ (-)0.14.0[39]
A eff (-)12.537.5
β 0 (-)1.154.60
β 1 (-)−3.6−1.2
S m a x ref (kPa)90140[61]
Land model (passive)
a p (kPa)1100
Rice model
[ C a 2 + ] res ( μ M)0.0750.225
[ C a 2 + ] max ( μ M)0.30.9
τ CR (ms)64.5192.5
τ CD (ms)64192

Appendix C. Biomarker Definitions

The definitions of the biomarkers used in this study are illustrated in Figure A1.
Figure A1. Biomarker definition and critical time points. (a) Pressure biomarkers, (b) volume biomarkers, and (c) active stress biomarkers.
Figure A1. Biomarker definition and critical time points. (a) Pressure biomarkers, (b) volume biomarkers, and (c) active stress biomarkers.
Mathematics 10 00823 g0a1
Based on the cavity pressure trace, the pressure biomarkers are defined as:
P max = max P ( t ) , d P d t | max = max t < t 2 d P d t , d P d t | min = min t > t 2 d P d t , P T D 90 = t 3 t 1 .
Here, P max is the maximum, d P d t | max and d P d t | min are the maximum and minimum rates, and P T D 90 is the transient duration at 90% decay from the maximum, measured between the time points at 10% of the amplitude.
Based on the cavity volume trace, the volume biomarkers are defined as:
S V = max V ( t ) min V ( t ) , τ VD = 1 s VD , τ VR = 1 s VR .
Here, S V is the stroke volume, τ VD is the decay time constant during the ejection phase, and τ VR is the rise time constants for the filling phase until the beginning of the atrial kick, which is defined to be the time point at 80% volume recovery, in line with [73]. The slope s VD is computed by an ordinary linear regression of
ln V ( t ) min V ( t ) S V
against t from t 4 to t 5 and the slope s VR is computed by an ordinary linear regression of
ln 1 V ( t ) min V ( t ) V ( t 7 ) min V ( t )
against t from t 6 to t 7 .
Based on the cellular active stress trace, the active stress biomarkers are defined as:
S a max = max S a ( t ) , d S a d t | max = max t < t 10 d S a d t , S T D 30 = t 8 t 12 , S T D 50 = t 9 t 13 , S T D 90 = t 10 t 14 .
Here, S a max is the maximum, d S a d t | max is the maximum rate, and S T D 30 , S T D 50 and S T D 90 are the transient durations at 30%, 50%, and 90% decay from the maximum, measured between the time points at 70%, 50%, and 10% of the amplitude.

Appendix D. Penalty Formulation of the Constrained Minimisation Problem in the Second Step of the Active Mechanics Personalisation Approach

The constrained minimisation problem Equation (19) was solved as a series of unconstrained problems by application of the penalty method:
min p HFM , a p , p CEM j = 1 n HFM w S a j B S a tar , j B S a sim , j ( p HFM , a p , p CEM ) B S a tar , j 2 + w λ B λ tar B λ sim ( p HFM , a p , p CEM ) B λ tar 2 + pen .
where the penalty term reads as follows:
pen = S a res 2 + 0 T [ min ( 0 , S a ( t ) ) ] 2 d t + t S a max T max 0 , d S a ( t ) d t 2 d t + [ min ( 0 , C T D 50 120 ) ] 2 + max ( 0 , C T D 50 120 ) ] 2 + [ min ( 0 , C T D 90 220 ) ] 2 + max ( 0 , C T D 90 785 ) ] 2 .
Since the differential evolution method was used for solving, the parameter constraints were enforced by limiting the populations of candidate solutions to the admissible ranges.

Appendix E. Estimated Parameter Values of the Personalised Models

The values of the Tanh model parameters and the values of the VFM parameters that were estimated in the first step are given in Table A6. The values of the Land model parameters and the Rice model parameters that were estimated in the second step are given in Table A7.
Table A6. Estimated parameter values of the personalised Tanh model and the personalised valve flow models for each patient case.
Table A6. Estimated parameter values of the personalised Tanh model and the personalised valve flow models for each patient case.
Patient CaseSmaxref τ SR ref τ SD TCRRAVfRMVf
(kPa)(ms)(ms)(ms)(mmHg mL s 1 )(mmHg mL s 1 )
01-CoA108.420.536.5376.40.00680.0444
02-CoA103.230.148.8529.30.01250.0794
03-CoA182.37.077.5514.80.00280.0562
04-CoA220.40.8935.0330.00.00010.0646
05-CoA145.212.539.2352.40.00400.0472
06-CoA172.16.140.3376.00.00260.0324
07-CoA146.826.164.9475.50.00690.0593
Table A7. Estimated parameter values of the personalised Land model and the personalised Rice model for each patient case.
Table A7. Estimated parameter values of the personalised Land model and the personalised Rice model for each patient case.
Patient CasenTRPN β 1 nTm[Ca2+]50refkUW[Ca2+]res[Ca2+]max τ CR τ CD
(-)(-)(-)( μ M)(ms 1 )( μ M)( μ M)(ms)(ms)
01-CoA3.88−1.203.010.6150.1470.0760.895105.6120.0
02-CoA3.22−1.215.470.5010.0460.0760.900164.4149.0
03-CoA4.99−1.206.020.5000.2500.0750.900162.9152.8
04-CoA5.00−1.205.680.5000.0310.0750.899140.467.6
05-CoA4.82−1.243.000.5050.0310.0820.898103.693.4
06-CoA4.99−1.203.480.5110.0290.0760.899113.8114.5
07-CoA4.98−1.204.120.5000.0970.0760.899139.4142.6

Appendix F. Convergence Behaviour for Patient Case 04-CoA in the First Step

Figure A2 illustrates the convergence behaviour in the first step of the personalisation approach for patient case 04-CoA. It shows that parameter values and cost both converge, however, after the ninth iteration, the required convergence threshold is still not reached. Ultimately, the updated set of parameter values after the ninth iteration caused the simulation to fail.
Figure A2. Convergence behaviour of the parameter value and the cost. For each iteration, the parameter values relative to their respective maxima (left y-axis; solid coloured lines) and the cost (right y-axis; solid black line) are plotted. In addition, the convergence threshold (0.1, dotted black line) is given.
Figure A2. Convergence behaviour of the parameter value and the cost. For each iteration, the parameter values relative to their respective maxima (left y-axis; solid coloured lines) and the cost (right y-axis; solid black line) are plotted. In addition, the convergence threshold (0.1, dotted black line) is given.
Mathematics 10 00823 g0a2

References

  1. Cardiovascular Diseases (CVDs). 2021. Available online: https://www.who.int/news-room/fact-sheets/detail/cardiovascular-diseases-(cvds) (accessed on 18 November 2021).
  2. Niederer, S.A.; Lumens, J.; Trayanova, N.A. Computational models in cardiology. Nat. Rev. Cardiol. 2018, 16, 100–111. [Google Scholar] [CrossRef] [PubMed]
  3. Niederer, S.A.; Campbell, K.S.; Campbell, S.G. A short history of the development of mathematical models of cardiac mechanics. J. Mol. Cell. Cardiol. 2019, 127, 11–19. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Trayanova, N. From genetics to smart watches: Developments in precision cardiology. Nat. Rev. Cardiol. 2018, 16, 72–73. [Google Scholar] [CrossRef]
  5. Corral-Acero, J.; Margara, F.; Marciniak, M.; Rodero, C.; Loncaric, F.; Feng, Y.; Gilbert, A.; Fernandes, J.F.; Bukhari, H.A.; Wajdan, A.; et al. The ‘Digital Twin’ to enable the vision of precision cardiology. Eur. Heart J. 2020, 41, 4556–4564. [Google Scholar] [CrossRef] [PubMed]
  6. Baillargeon, B.; Rebelo, N.; Fox, D.D.; Taylor, R.L.; Kuhl, E. The Living Heart Project: A robust and integrative simulator for human heart function. Eur. J. Mech.-A/Solids 2014, 48, 38–47. [Google Scholar] [CrossRef] [Green Version]
  7. Augustin, C.M.; Neic, A.; Liebmann, M.; Prassl, A.J.; Niederer, S.A.; Haase, G.; Plank, G. Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation. J. Comput. Phys. 2016, 305, 622–646. [Google Scholar] [CrossRef] [Green Version]
  8. Strocchi, M.; Gsell, M.A.; Augustin, C.M.; Razeghi, O.; Roney, C.H.; Prassl, A.J.; Vigmond, E.J.; Behar, J.M.; Gould, J.S.; Rinaldi, C.A.; et al. Simulating ventricular systolic motion in a four-chamber heart model with spatially varying robin boundary conditions to model the effect of the pericardium. J. Biomech. 2020, 101, 109645. [Google Scholar] [CrossRef]
  9. Regazzoni, F.; Salvador, M.; Africa, P.C.; Fedele, M.; Dede’, L.; Quarteroni, A. A cardiac electromechanics model coupled with a lumped parameters model for closed-loop blood circulation. Part II: Numerical approximation. arXiv 2020, arXiv:2011.15051. [Google Scholar]
  10. Levrero-Florencio, F.; Margara, F.; Zacur, E.; Bueno-Orovio, A.; Wang, Z.; Santiago, A.; Aguado-Sierra, J.; Houzeaux, G.; Grau, V.; Kay, D.; et al. Sensitivity analysis of a strongly-coupled human-based electromechanical cardiac model: Effect of mechanical parameters on physiologically relevant biomarkers. Comput. Methods Appl. Mech. Eng. 2020, 361, 112762. [Google Scholar] [CrossRef]
  11. Augustin, C.M.; Gsell, M.A.; Karabelas, E.; Willemen, E.; Prinzen, F.W.; Lumens, J.; Vigmond, E.J.; Plank, G. A computationally efficient physiologically comprehensive 3D–0D closed-loop model of the heart and circulation. Comput. Methods Appl. Mech. Eng. 2021, 386, 114092. [Google Scholar] [CrossRef]
  12. Gerach, T.; Schuler, S.; Fröhlich, J.; Lindner, L.; Kovacheva, E.; Moss, R.; Wülfers, E.M.; Seemann, G.; Wieners, C.; Loewe, A. Electro-Mechanical Whole-Heart Digital Twins: A Fully Coupled Multi-Physics Approach. Mathematics 2021, 9, 1247. [Google Scholar] [CrossRef]
  13. Mirams, G.R.; Pathmanathan, P.; Gray, R.A.; Challenor, P.; Clayton, R.H. Uncertainty and variability in computational and mathematical models of cardiac physiology. J. Physiol. 2016, 594, 6833–6847. [Google Scholar] [CrossRef] [PubMed]
  14. Niederer, S.A.; Sacks, M.S.; Girolami, M.; Willcox, K. Scaling digital twins from the artisanal to the industrial. Nat. Comput. Sci. 2021, 1, 313–320. [Google Scholar] [CrossRef]
  15. Camps, J.; Lawson, B.; Drovandi, C.; Minchole, A.; Wang, Z.J.; Grau, V.; Burrage, K.; Rodriguez, B. Inference of ventricular activation properties from non-invasive electrocardiography. Med. Image Anal. 2021, 73, 102143. [Google Scholar] [CrossRef]
  16. Gillette, K.; Gsell, M.A.; Prassl, A.J.; Karabelas, E.; Reiter, U.; Reiter, G.; Grandits, T.; Payer, C.; Štern, D.; Urschler, M.; et al. A Framework for the generation of digital twins of cardiac electrophysiology from clinical 12-leads ECGs. Med. Image Anal. 2021, 71, 102080. [Google Scholar] [CrossRef] [PubMed]
  17. Pezzuto, S.; Prinzen, F.W.; Potse, M.; Maffessanti, F.; Regoli, F.; Caputo, M.L.; Conte, G.; Krause, R.; Auricchio, A. Reconstruction of three-dimensional biventricular activation based on the 12-lead electrocardiogram via patient-specific modelling. Europace 2021, 23, 640–647. [Google Scholar] [CrossRef]
  18. Marx, L.; Gsell, M.; Rund, A.; Caforio, F.; Prassl, A.J.; Toth-Gayor, G.; Kuehne, T.; Augustin, C.; Plank, G. Personalization of electro-mechanical models of the pressure-overloaded left ventricle: Fitting of Windkessel-type afterload models. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2020, 378, 20190342. [Google Scholar] [CrossRef]
  19. Cai, L.; Ren, L.; Wang, Y.; Xie, W.; Zhu, G.; Gao, H. Surrogate models based on machine learning methods for parameter estimation of left ventricular myocardium. R. Soc. Open Sci. 2021, 8, 201121. [Google Scholar] [CrossRef]
  20. Marx, L.; Niestrawska, J.A.; Gsell, M.A.F.; Caforio, F.; Plank, G.; Augustin, C.M. Efficient identification of myocardial material parameters and the stress-free reference configuration for patient-specific human heart models. arXiv 2021, arXiv:q-bio.TO/2101.04411. [Google Scholar]
  21. Ambrosi, D.; Pezzuto, S. Active Stress vs. Active Strain in Mechanobiology: Constitutive Issues. J. Elast. 2012, 107, 199–212. [Google Scholar] [CrossRef]
  22. Bers, D. Excitation-Contraction Coupling and Cardiac Contractile Force, 2nd ed.; Springer: Dordrecht, The Netherlands, 2001. [Google Scholar] [CrossRef]
  23. Kayvanpour, E.; Mansi, T.; Sedaghat-Hamedani, F.; Amr, A.; Neumann, D.; Georgescu, B.; Seegerer, P.; Kamen, A.; Haas, J.; Frese, K.S.; et al. Towards Personalized Cardiology: Multi-Scale Modeling of the Failing Heart. PLoS ONE 2015, 10, e0134869. [Google Scholar] [CrossRef] [PubMed]
  24. Asner, L.; Hadjicharalambous, M.; Chabiniok, R.; Peresutti, D.; Sammut, E.; Wong, J.; Carr-White, G.; Chowienczyk, P.; Lee, J.; King, A.; et al. Estimation of passive and active properties in the human heart using 3D tagged MRI. Biomech. Model. Mechanobiol. 2016, 15, 1121–1139. [Google Scholar] [CrossRef] [PubMed]
  25. Finsberg, H.; Xi, C.; Tan, J.L.; Zhong, L.; Genet, M.; Sundnes, J.; Lee, L.C.; Wall, S.T. Efficient estimation of personalized biventricular mechanical function employing gradient-based optimization. Int. J. Numer. Methods Biomed. Eng. 2018, 34, e2982. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Longobardi, S.; Sher, A.; Niederer, S.A. In silico identification of potential calcium dynamics and sarcomere targets for recovering left ventricular function in rat heart failure with preserved ejection fraction. PLoS Comput. Biol. 2021, 17, e1009646. [Google Scholar] [CrossRef]
  27. Longobardi, S.; Lewalle, A.; Coveney, S.; Sjaastad, I.; Espe, E.K.S.; Louch, W.E.; Musante, C.J.; Sher, A.; Niederer, S.A. Predicting left ventricular contractile function via Gaussian process emulation in aortic-banded rats. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2020, 378, 20190334. [Google Scholar] [CrossRef]
  28. Fernandes, J.F.; Goubergrits, L.; Brüning, J.; Hellmeier, F.; Nordmeyer, S.; da Silva, T.F.; Schubert, S.; Berger, F.; Kuehne, T.; Kelm, M.; et al. Beyond pressure gradients: The effects of intervention on heart power in aortic coarctation. PLoS ONE 2017, 12, e0168487. [Google Scholar] [CrossRef] [Green Version]
  29. Bols, J.; Degroote, J.; Trachet, B.; Verhegghe, B.; Segers, P.; Vierendeels, J. A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels. J. Comput. Appl. Math. 2013, 246, 10–17. [Google Scholar] [CrossRef]
  30. Bayer, J.; Blake, R.; Plank, G.; Trayanova, N. A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 2012, 40, 2243–2254. [Google Scholar] [CrossRef] [Green Version]
  31. Prassl, A.; Kickinger, F.; Ahammer, H.; Grau, V.; Schneider, J.; Hofer, E.; Vigmond, E.; Trayanova, N.; Plank, G. Automatically Generated, Anatomically Accurate Meshes for Cardiac Electrophysiology Problems. IEEE Trans. Biomed. Eng. 2009, 56, 1318–1330. [Google Scholar] [CrossRef] [Green Version]
  32. Crozier, A.; Augustin, C.; Neic, A.; Prassl, A.; Holler, M.; Fastl, T.; Hennemuth, A.; Bredies, K.; Kuehne, T.; Bishop, M.; et al. Image-based personalization of cardiac anatomy for coupled electromechanical modeling. Ann. Biomed. Eng. 2016, 44, 58–70. [Google Scholar] [CrossRef] [Green Version]
  33. Flory, P.J. Thermodynamic relations for high elastic materials. Trans. Faraday Soc. 1961, 57, 829–838. [Google Scholar] [CrossRef]
  34. Augustin, C.M.; Crozier, A.; Neic, A.; Prassl, A.J.; Karabelas, E.; Ferreira da Silva, T.; Fernandes, J.F.; Campos, F.; Kuehne, T.; Plank, G. Patient-specific modeling of left ventricular electromechanics as a driver for haemodynamic analysis. EP Eur. 2016, 18, iv121–iv129. [Google Scholar] [CrossRef] [PubMed]
  35. Westerhof, N.; Lankhaar, J.; Westerhof, B. The arterial Windkessel. Med. Biol. Eng. Comput. 2009, 47, 131–141. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Guccione, J.; McCulloch, A.; Waldman, L. Passive material properties of intact ventricular myocardium determined from a cylindrical model. J. Biomech. Eng. 1991, 113, 42–55. [Google Scholar] [CrossRef] [PubMed]
  37. Kerckhoffs, R.; Bovendeert, P.; Prinzen, F.; Smits, K.; Arts, T. Intra- and interventricular asynchrony of electromechanics in the ventricularly paced heart. J. Eng. Math. 2003, 47, 201–216. [Google Scholar] [CrossRef]
  38. Niederer, S.A.; Plank, G.; Chinchapatnam, P.; Ginks, M.; Lamata, P.; Rhode, K.S.; Rinaldi, C.A.; Razavi, R.; Smith, N.P. Length-dependent tension in the failing heart and the efficacy of cardiac resynchronization therapy. Cardiovasc. Res. 2011, 89, 336–343. [Google Scholar] [CrossRef] [Green Version]
  39. Land, S.; Park-Holohan, S.J.; Smith, N.P.; dos Remedios, C.G.; Kentish, J.C.; Niederer, S.A. A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. J. Mol. Cell. Cardiol. 2017, 106, 68–83. [Google Scholar] [CrossRef] [Green Version]
  40. Neic, A.; Campos, F.; Prassl, A.; Niederer, S.; Bishop, M.; Vigmond, E.; Plank, G. Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model. J. Comput. Phys. 2017, 346, 191–211. [Google Scholar] [CrossRef]
  41. Luo, C.; Rudy, Y. A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circ. Res. 1991, 68, 1501–1526. [Google Scholar] [CrossRef] [Green Version]
  42. Rice, J.J.; Wang, F.; Bers, D.M.; De Tombe, P.P. Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J. 2008, 95, 2368–2390. [Google Scholar] [CrossRef] [Green Version]
  43. Coppini, R.; Ferrantini, C.; Yao, L.; Fan, P.; Del Lungo, M.; Stillitano, F.; Sartiani, L.; Tosi, B.; Suffredini, S.; Tesi, C.; et al. Late sodium current inhibition reverses electromechanical dysfunction in human hypertrophic cardiomyopathy. Circulation 2013, 127, 575–584. [Google Scholar] [CrossRef] [PubMed]
  44. Sutanto, H.; Lyon, A.; Lumens, J.; Schotten, U.; Dobrev, D.; Heijman, J. Cardiomyocyte calcium handling in health and disease: Insights from in vitro and in silico studies. Prog. Biophys. Mol. Biol. 2020, 157, 54–75. [Google Scholar] [CrossRef]
  45. Quinn, T.; Kohl, P. Cardiac mechano-electric coupling: Acute effects of mechanical stimulation on heart rate and rhythm. Physiol. Rev. 2021, 101, 37–92. [Google Scholar] [CrossRef] [PubMed]
  46. Regazzoni, F.; Quarteroni, A. An oscillation-free fully staggered algorithm for velocity-dependent active models of cardiac mechanics. Comput. Methods Appl. Mech. Eng. 2021, 373, 113506. [Google Scholar] [CrossRef]
  47. Vigmond, E.J.; dos Santos, R.; Prassl, A.J.; Deo, M.; Plank, G.; Weber dos Santos, R.; Prassl, A.J.; Deo, M.; Plank, G. Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 2008, 96, 3–18. [Google Scholar] [CrossRef] [Green Version]
  48. Plank, G.; Loewe, A.; Neic, A.; Augustin, C.; Huang, Y.L.; Gsell, M.A.; Karabelas, E.; Nothstein, M.; Prassl, A.J.; Sánchez, J.; et al. The openCARP simulation environment for cardiac electrophysiology. Comput. Methods Programs Biomed. 2021, 208, 106223. [Google Scholar] [CrossRef]
  49. Neic, A.; Liebmann, M.; Hoetzl, E.; Mitchell, L.; Vigmond, E.; Haase, G.; Plank, G. Accelerating cardiac bidomain simulations using graphics processing units. IEEE Trans. Biomed. Eng. 2012, 59, 2281–2290. [Google Scholar] [CrossRef] [Green Version]
  50. Niederer, S.; Kerfoot, E.; Benson, A.; Bernabeau, M.; Bernus, O.; Bradley, C.; Cherry, E.; Clayton, R.; Fenton, F.; Garny, A.; et al. Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2011, 369, 4331–4351. [Google Scholar] [CrossRef] [Green Version]
  51. Land, S.; Gurev, V.; Arens, S.; Augustin, C.; Baron, L.; Blake, R.; Bradley, C.; Castro, S.; Crozier, A.; Favino, M.; et al. Verification of cardiac mechanics software: Benchmark problems and solutions for testing active and passive material behaviour. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2015, 471, 20150641. [Google Scholar] [CrossRef] [Green Version]
  52. Passini, E.; Trovato, C.; Morissette, P.; Sannajust, F.; Bueno-Orovio, A.; Rodriguez, B. Drug-induced shortening of the electromechanical window is an effective biomarker for in silico prediction of clinical risk of arrhythmias. Br. J. Pharmacol. 2019, 176, 3819–3833. [Google Scholar] [CrossRef]
  53. Storn, R.; Price, K. Differential Evolution—A Simple and Efficient Heuristic for global Optimization over Continuous Spaces. J. Glob. Optim. 1997, 11, 341–359. [Google Scholar] [CrossRef]
  54. Vesterstrom, J.; Thomsen, R. A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems. In Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753), Portland, OR, USA, 19–23 June 2004. [Google Scholar] [CrossRef]
  55. Rocca, P.; Oliveri, G.; Massa, A. Differential evolution as applied to electromagnetics. IEEE Antennas Propag. Mag. 2011, 53, 38–49. [Google Scholar] [CrossRef]
  56. Newville, M.; Stensitzki, T.; Allen, D.; Ingargiola, A. LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python. Astrophys. Source Code Libr. 2016. [Google Scholar] [CrossRef]
  57. Sobol, I.M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simul. 2001, 55, 271–280. [Google Scholar] [CrossRef]
  58. Saltelli, A.; Ratto, M.; Andres, T.; Campolongo, F.; Cariboni, J.; Gatelli, D.; Saisana, M.; Tarantola, S. Global Sensitivity Analysis. The Primer; John Wiley and Sons: Hoboken, NJ, USA, 2008; pp. 1–292. [Google Scholar] [CrossRef] [Green Version]
  59. Saltelli, A. Making best use of model evaluations to compute sensitivity indices. Comput. Phys. Commun. 2002, 145, 280–297. [Google Scholar] [CrossRef]
  60. Herman, J.; Usher, W. SALib: An open-source Python library for Sensitivity Analysis. J. Open Source Softw. 2017, 2, 97. [Google Scholar] [CrossRef]
  61. Tøndel, K.; Land, S.; Niederer, S.A.; Smith, N.P. Quantifying inter-species differences in contractile function through biophysical modelling. J. Physiol. 2015, 593, 1083–1111. [Google Scholar] [CrossRef] [Green Version]
  62. Sermesant, M.; Delingette, H.; Ayache, N. An electromechanical model of the heart for image analysis and simulation. IEEE Trans. Med. Imaging 2006, 25, 612–625. [Google Scholar] [CrossRef]
  63. Russell, K.; Eriksen, M.; Aaberge, L.; Wilhelmsen, N.; Skulstad, H.; Remme, E.W.; Haugaa, K.H.; Opdahl, A.; Fjeld, J.G.; Gjesdal, O.; et al. A novel clinical method for quantification of regional left ventricular pressure–strain loop area: A non-invasive index of myocardial work. Eur. Heart J. 2012, 33, 724–733. [Google Scholar] [CrossRef] [Green Version]
  64. Land, S.; Niederer, S.A.; Aronsen, J.M.; Espe, E.K.S.; Zhang, L.; Louch, W.E.; Sjaastad, I.; Sejersted, O.M.; Smith, N.P. An analysis of deformation-dependent electromechanical coupling in the mouse heart. J. Physiol. 2012, 590, 4553–4569. [Google Scholar] [CrossRef]
  65. Nagueh, S.F.; Middleton, K.J.; Kopelen, H.A.; Zoghbi, W.A.; Quiñones, M.A. Doppler Tissue Imaging: A Noninvasive Technique for Evaluation of Left Ventricular Relaxation and Estimation of Filling Pressures. J. Am. Coll. Cardiol. 1997, 30, 1527–1533. [Google Scholar] [CrossRef]
  66. Nagueh, S.F.; Appleton, C.P.; Gillebert, T.C.; Marino, P.N.; Oh, J.K.; Smiseth, O.A.; Waggoner, A.D.; Flachskampf, F.A.; Pellikka, P.A.; Evangelisa, A. Recommendations for the Evaluation of Left Ventricular Diastolic Function by Echocardiography. Eur. J. Echocardiogr. 2008, 10, 165–193. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Brett, S.E.; Guilcher, A.; Clapp, B.; Chowienczyk, P. Estimating central systolic blood pressure during oscillometric determination of blood pressure. Blood Press. Monit. 2012, 17, 132–136. [Google Scholar] [CrossRef] [PubMed]
  68. Alber, M.; Tepole, A.B.; Cannon, W.R.; De, S.; Dura-Bernal, S.; Garikipati, K.; Karniadakis, G.; Lytton, W.W.; Perdikaris, P.; Petzold, L.; et al. Integrating machine learning and multiscale modeling—Perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences. npj Digit. Med. 2019, 2, 115. [Google Scholar] [CrossRef]
  69. Lin, D.H.S.; Yin, F.C.P. A Multiaxial Constitutive Law for Mammalian Left Ventricular Myocardium in Steady-State Barium Contracture or Tetanus. J. Biomech. Eng. 1998, 120, 504–517. [Google Scholar] [CrossRef]
  70. Walker, J.C.; Ratcliffe, M.B.; Zhang, P.; Wallace, A.W.; Fata, B.; Hsu, E.W.; Saloner, D.; Guccione, J.M. MRI-based finite-element analysis of left ventricular aneurysm. Am. J.-Physiol.-Heart Circ. Physiol. 2005, 289, H692–H700. [Google Scholar] [CrossRef] [Green Version]
  71. Genet, M.; Lee, L.C.; Nguyen, R.; Haraldsson, H.; Acevedo-Bolton, G.; Zhang, Z.; Ge, L.; Ordovas, K.; Kozerke, S.; Guccione, J.M. Distribution of normal human left ventricular myofiber stress at end diastole and end systole: A target for in silico design of heart failure treatments. J. Appl. Physiol. 2014, 117, 142–152. [Google Scholar] [CrossRef]
  72. Tomek, J.; Bueno-Orovio, A.; Passini, E.; Zhou, X.; Minchole, A.; Britton, O.; Bartolucci, C.; Severi, S.; Shrier, A.; Virag, L.; et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. eLife 2019, 8, e48890. [Google Scholar] [CrossRef]
  73. Namana, V.; Gupta, S.; Sabharwal, N.; Hollander, G. Clinical significance of atrial kick. QJM Int. J. Med. 2018, 111, 569–570. [Google Scholar] [CrossRef] [Green Version]
Figure 1. (Left), seven patient-specific anatomical models of the left ventricle (LV) and the aortic root from patients treated for aortic coarctation (CoA); (right), mechanical and afterload boundary conditions where Neumann-type pressure boundary conditions are illustrated in blue and Robin-type boundary conditions are illustrated in green (uni-directional springs that mimic the effect of the pericardium) and yellow (omni-directional springs). The mechanical model is coupled to a 0D three-element Windkessel model (WK3) that accounts for afterload conditions. Here, R 1 , R 2 are the characteristic and peripheral resistances, respectively, and C is the arterial compliance. Pressure is denoted by P and the relationship between pressure and flow q = d V d t across the aortic valve (AV) is represented by a 0D diode model, where R AV is the respective resistance. The diode model of the flow across the mitral valve and the uni-directional springs that are applied to the septum are not shown. See [18] for more details.
Figure 1. (Left), seven patient-specific anatomical models of the left ventricle (LV) and the aortic root from patients treated for aortic coarctation (CoA); (right), mechanical and afterload boundary conditions where Neumann-type pressure boundary conditions are illustrated in blue and Robin-type boundary conditions are illustrated in green (uni-directional springs that mimic the effect of the pericardium) and yellow (omni-directional springs). The mechanical model is coupled to a 0D three-element Windkessel model (WK3) that accounts for afterload conditions. Here, R 1 , R 2 are the characteristic and peripheral resistances, respectively, and C is the arterial compliance. Pressure is denoted by P and the relationship between pressure and flow q = d V d t across the aortic valve (AV) is represented by a 0D diode model, where R AV is the respective resistance. The diode model of the flow across the mitral valve and the uni-directional springs that are applied to the septum are not shown. See [18] for more details.
Mathematics 10 00823 g001
Figure 2. Conceptual illustration of the two-step multi-fidelity approach for personalising biophysically detailed active mechanics models. In the first step, the active mechanical behaviour in the given 3D EM model is described by a low-fidelity model (LFM). This is personalised at the organ scale by minimising the difference between simulated (sim) and clinical (clin) pressure biomarker values ( B P ) that are obtained from an available clinical cavity pressure trace. Median traces of nodal cellular active stress ( S a ), intracellular calcium concentration ( [ C a 2 + ] i ), and fibre stretch ( λ ) are then obtained from the simulation that was produced by the personalised 3D EM models. These are utilised to personalise the high-fidelity model (HFM) at the cell scale in the second step. To this end, the HFM model is integrated in an 0D EM model that simulates the stretch of the cardiomyocyte during the cardiac cycle based on the equilibrium of active and passive stress. The median trace of nodal [ C a 2 + ] i is used as input but [ C a 2 + ] i can also be generated from a coupled model of [ C a 2 + ] i evolution that is integrated in the chosen cellular EP model. Personalisation is done by minimising the difference between simulated (sim) and target (tar) biomarker values that are obtained from the median traces of nodal cellular active stress ( B S a ) and fibre stretch ( B λ ). The parameters of the LFM and the HFM are denoted by p LFM and p HFM , respectively. Please note that the fibre stretch in the relaxed state is 1 per definition (Section 2.3.1 and Appendix A.3).
Figure 2. Conceptual illustration of the two-step multi-fidelity approach for personalising biophysically detailed active mechanics models. In the first step, the active mechanical behaviour in the given 3D EM model is described by a low-fidelity model (LFM). This is personalised at the organ scale by minimising the difference between simulated (sim) and clinical (clin) pressure biomarker values ( B P ) that are obtained from an available clinical cavity pressure trace. Median traces of nodal cellular active stress ( S a ), intracellular calcium concentration ( [ C a 2 + ] i ), and fibre stretch ( λ ) are then obtained from the simulation that was produced by the personalised 3D EM models. These are utilised to personalise the high-fidelity model (HFM) at the cell scale in the second step. To this end, the HFM model is integrated in an 0D EM model that simulates the stretch of the cardiomyocyte during the cardiac cycle based on the equilibrium of active and passive stress. The median trace of nodal [ C a 2 + ] i is used as input but [ C a 2 + ] i can also be generated from a coupled model of [ C a 2 + ] i evolution that is integrated in the chosen cellular EP model. Personalisation is done by minimising the difference between simulated (sim) and target (tar) biomarker values that are obtained from the median traces of nodal cellular active stress ( B S a ) and fibre stretch ( B λ ). The parameters of the LFM and the HFM are denoted by p LFM and p HFM , respectively. Please note that the fibre stretch in the relaxed state is 1 per definition (Section 2.3.1 and Appendix A.3).
Mathematics 10 00823 g002
Figure 3. Global sensitivity of active stress biomarkers to parameters of the Land and the Rice model in the 0D EM model. The Sobol global sensitivity analysis was performed and first-order ( S 1 ) and total-order sensitivity indices ( S T ) are shown.
Figure 3. Global sensitivity of active stress biomarkers to parameters of the Land and the Rice model in the 0D EM model. The Sobol global sensitivity analysis was performed and first-order ( S 1 ) and total-order sensitivity indices ( S T ) are shown.
Mathematics 10 00823 g003
Figure 4. Comparison of simulated and clinical left ventricular pressure and volume after personalisation. The blue/red lines represent the simulation data Sim1/Sim2 that were produced by the 3D model of human left ventricular EM in the first/second step of the active mechanics personalisation approach. The black dots represent clinical data (Clin). (a) Pressure–volume loops, (b) pressure traces, and (c) volume traces are shown.
Figure 4. Comparison of simulated and clinical left ventricular pressure and volume after personalisation. The blue/red lines represent the simulation data Sim1/Sim2 that were produced by the 3D model of human left ventricular EM in the first/second step of the active mechanics personalisation approach. The black dots represent clinical data (Clin). (a) Pressure–volume loops, (b) pressure traces, and (c) volume traces are shown.
Mathematics 10 00823 g004
Figure 5. Comparison of simulated cellular active stress and intracellular calcium concentration. The blue lines represent the median nodal traces obtained from the simulations that were produced by the 3D model of human left ventricular EM in the first step (Sim1). The red lines represent the traces that were produced by the 0D model of cellular EM in the second step of the active mechanics personalisation approach (Sim2). (a) Cellular active stress traces (solid line represents a target for the personalisation), (b) intracellular calcium concentration traces (dashed line does not represent a target for the personalisation).
Figure 5. Comparison of simulated cellular active stress and intracellular calcium concentration. The blue lines represent the median nodal traces obtained from the simulations that were produced by the 3D model of human left ventricular EM in the first step (Sim1). The red lines represent the traces that were produced by the 0D model of cellular EM in the second step of the active mechanics personalisation approach (Sim2). (a) Cellular active stress traces (solid line represents a target for the personalisation), (b) intracellular calcium concentration traces (dashed line does not represent a target for the personalisation).
Mathematics 10 00823 g005
Figure 6. Effect of clinical data uncertainty on simulated left ventricular pressure and volume in the two steps of the active mechanics personalisation approach. The blue/red lines represent the simulation data (Sim1, Sim2) that were produced by the 3D model of human left ventricular EM in the first/second step of the active mechanics personalisation approach. The active mechanics personalisation was performed based on ten samples of clinical biomarkers (±10% around the measured value) and the resulting pressures and volumes (light colours) are compared to those that resulted from the personalisation based on the measured values (bold colours). (a) Pressure–volume loops, (b) pressure traces, (c) volume traces.
Figure 6. Effect of clinical data uncertainty on simulated left ventricular pressure and volume in the two steps of the active mechanics personalisation approach. The blue/red lines represent the simulation data (Sim1, Sim2) that were produced by the 3D model of human left ventricular EM in the first/second step of the active mechanics personalisation approach. The active mechanics personalisation was performed based on ten samples of clinical biomarkers (±10% around the measured value) and the resulting pressures and volumes (light colours) are compared to those that resulted from the personalisation based on the measured values (bold colours). (a) Pressure–volume loops, (b) pressure traces, (c) volume traces.
Mathematics 10 00823 g006
Table 1. Mesh information, CPU times, and iteration numbers for each patient case. Given is the number of nodes and elements of the patient-specific FE mesh, the cardiac cycle length, and the times per iteration and number of iterations for personalising the low-fidelity model (LFM) in the first step (3D EM model) and the high fidelity model (HFM) in the second step (0D EM model). 3D simulations in the first step were run on ARCHER2 using 128 cores and 0D simulations in the second step were run on a desktop computer using 1 core.
Table 1. Mesh information, CPU times, and iteration numbers for each patient case. Given is the number of nodes and elements of the patient-specific FE mesh, the cardiac cycle length, and the times per iteration and number of iterations for personalising the low-fidelity model (LFM) in the first step (3D EM model) and the high fidelity model (HFM) in the second step (0D EM model). 3D simulations in the first step were run on ARCHER2 using 128 cores and 0D simulations in the second step were run on a desktop computer using 1 core.
Patient Case# Nodes# ElementsCycle LengthTime/It LFMTime/It HFM# It LFM# It HFM
(-)(-)(ms)(s)(s)(-)(-)
01-CoA159,948806,4306594834.70.033213,211
02-CoA162,188835,51612317945.40.061316,061
03-CoA63,804301,1469172389.40.045417,261
04-CoA96,176487,1326312249.10.03198561
05-CoA126,981652,0126543434.80.032316,211
06-CoA165,508853,7176975407.90.036414,111
07-CoA82,212394,6908522591.80.042320,766
Table 2. Simulated and clinical LV pressure and volume biomarker values for each patient case. Clinical data are given in the first panel and results of the first and the second step of the active mechanics personalisation approach are given in the second and third panel. The goodness of fit is measured as relative difference between simulated and clinical value and given in brackets.
Table 2. Simulated and clinical LV pressure and volume biomarker values for each patient case. Clinical data are given in the first panel and results of the first and the second step of the active mechanics personalisation approach are given in the second and third panel. The goodness of fit is measured as relative difference between simulated and clinical value and given in brackets.
Patient CasePmax d P d t | max d P d t | min PTD90SV
(mmHg)(mmHg ms 1 )(mmHg ms 1 )(ms)(mL)
1192.7−2.9309100
01-CoA115 (3.0%)1.9 (30.9%)−2.2 (26.7%)330 (6.7%)90 (10.4%)
125 (5.1%)2.8 (3.7%)−1.2 (59.4%)316 (2.1%)92.6 (7.4%)
1051.3−1.5452115
02-CoA105 (<0.1%)1.2 (5.3%)−1.5 (2.5%)460 (1.9%)123 (6.8%)
107 (1.5%)1.5 (18.5%)−0.9 (39.8%)446 (1.2%)126 (8.6%)
1211.9−1.245046
03-CoA114 (5.7%)1.7 (12.7%)−1.3 (5.4%)448 (0.3%)46 (0.1%)
112 (7.8%)2.5 (28.7%)−1.2 (1.2%)431 (4.1%)44 (4.0%)
1293.6−3.129068
04-CoA126 (2.7%)3.2 (12.0%)−3.1 (0.9%)290 (<0.1%)53 (20.9%)
126 (2.6%)3.8 (5.7%)−1.3 (55.8%)302 (4.1%)52 (23.6%)
1332.7−2.330992
05-CoA120 (10.1%)2.4 (13.8%)−2.5 (5.7%)309 (0.2%)81 (11.5%)
123 (7.4%)3.4 (22.9%)−1.1 (54.1%)310 (0.5%)83 (10.1%)
1523.2−2.7332100
06-CoA139 (8.7%)2.7 (15.2%)−2.8 (6.3%)332 (0.1%)87 (12.5%)
142 (6.4%)3.4 (6.8%)−1.2 (54.4%)345 (4.1%)83 (16.5%)
1101.2−1.241255
07-CoA 101 (7.7%)1.3 (5.2%)−1.2 (6.3%)409 (0.6%)60 (7.6%)
101 (7.5%)2.1 (71.2%)−1.0 (14.7%)394 (4.3%)59 (7.5%)
Table 3. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. Estimated model parameter values of the first step of the active mechanics personalisation approach are given. The original values are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 3. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. Estimated model parameter values of the first step of the active mechanics personalisation approach are given. The original values are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Smaxref τ SR ref τ SD TCRRAVfRMVf
(kPa)(ms)(ms)(ms)(mmHg mL s 1 )(mmHg mL s 1 )
Original103.230.148.8529.30.01250.0794
Mean103.430.748.3524.90.01120.0795
(7.8%)(8.0%)(4.6%)(5.0%)(10.8%)(5.0%)
SD9.42.62.730.30.00150.0048
(4.8%)(3.2%)(3.2%)(2.9%)(8.1%)(3.3%)
Min88.926.543.7478.90.00900.0728
(0.9%)(2.8%)(0.1%)(0.4%)(0.4%)(<0.1%)
Max118.833.651.9572.20.01380.0878
(15.1%)(13.9%)(10.4%)(9.5%)(23.0%)(10.6%)
Table 4. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. LV pressure and volume biomarker values of the first step of the active mechanics personalisation approach are given. The original values are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 4. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. LV pressure and volume biomarker values of the first step of the active mechanics personalisation approach are given. The original values are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Pmax d P d t | max d P d t | min PTD90SV
(mmHg)(mmHg ms 1 )(mmHg ms 1 )(ms)(mL)
Original1051.2−1.5460123
Mean1051.2−1.5457122
(3.6%)(6.5%)(4.7%)(5.4%)(4.3%)
SD40.10.1286
(2.2%)(3.4%)(2.5%)(2.9%)(2.8%)
Min991.1−1.6416112
(0.3%)(1.5%)(0.7%)(0.9%)(0.2%)
Max1121.3−1.4502132
(6.5%)(12.0%)(9.2%)(9.6%)(9.2%)
Table 5. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. Estimated model parameter values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 5. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. Estimated model parameter values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
nTRPN β 1 nTm[Ca2+]50refkUW[Ca2+]res[Ca2+]max τ CR τ CD
(-)(-)(-)( μ M)(ms 1 )( μ M)( μ M)(ms)(ms)
Original3.22−1.215.470.5010.0460.0760.900164.4149.0
Mean3.60−1.204.940.5070.0550.0750.900158.1154.0
(27.0%)(0.3%)(11.5%)(1.3%)(32.7%)(1.6%)(<0.1%)(5.4%)(4.8%)
SD0.94<0.010.760.0170.022<0.001<0.0017.06.2
(16.4%)(0.2%)(12.4%)(3.3%)(40.4%)(0.4%)(<0.1%)(2.1%)(2.4%)
Min2.42−1.213.010.5000.0360.0750.899147.7141.7
(9.1%)(0.0%)(1.8%)(0.1%)(1.1%)(0.8%)(<0.1%)(2.0%)(0.3%)
Max5.00−1.205.790.5570.1090.0760.900173.5162.2
(55.1%)(0.5%)(44.9%)(11.2%)(138.8%)(2.0%)(0.1%)(10.2%)(8.8%)
Table 6. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. LV pressure and volume biomarker values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 6. Results of the clinical data uncertainty robustness analysis for the patient case 02-CoA. LV pressure and volume biomarker values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of ten samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Pmax d P d t | max d P d t | min PTD90SV
(mmHg)[mmHg ms 1 ](mmHg ms 1 )(ms)(mL)
Original1071.5−0.9446126
Mean1081.6−1.0442124
(5.2%)(10.5%)(7.0%)(4.0%)(4.3%)
SD60.20.1206
(3.6%)(7.3%)(4.6%)(2.3%)(2.7%)
Min1001.3−1.1412115
(0.2%)(0.5%)(1.2%)(0.2%)(0.2%)
Max1191.8−0.8473135
(11.3%)(23.2%)(16.8%)(7.6%)(8.6%)
Table 7. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Estimated model parameters of the first step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 7. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Estimated model parameters of the first step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Smaxref τ SR ref τ SD TCRRAVfRMVf
[kPa][ms][ms][ms][mmHg mL s 1 ][mmHg mL s 1 ]
Original103.230.148.8529.30.01250.0794
Mean102.327.648.0517.60.01420.0817
(1.2%)(10.9%)(2.0%)(2.2%)(30.4%)(2.9%)
SD1.43.70.92.30.00360.0001
(1.1%)(11.5%)(1.5%)(0.4%)(28.5%)(0.1%)
Min100.420.746.9514.10.01010.0815
(0.2%)(1.1%)(0.1%)(1.7%)(0.5%)(2.7%)
Max104.031.249.2520.20.01920.0817
(2.7%)(33.0%)(3.9%)(2.9%)(70.6%)(3.0%)
Table 8. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Pressure and volume biomarker values of the first step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 8. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Pressure and volume biomarker values of the first step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Pmax d P d t | max d P d t | min PTD90SV
[mmHg][mmHg ms 1 ][mmHg ms 1 ][ms][mL]
Original1051.2−1.5460123
Mean1061.2−1.5451122
(1.0%)(3.4%)(3.1%)(1.9%)(1.4%)
SD106<0.1<0.111
(0.5%)(3.6%)(2.6%)(0.3%)(0.6%)
Min1041.2−1.6449121
(0.3%)(0.5%)(0.9%)(1.5%)(0.3%)
Max1071.3−1.5453123
(1.6%)(10.4%)(8.1%)(2.4%)(1.9%)
Table 9. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Estimated model parameters of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 9. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Estimated model parameters of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
nTRPN β 1 nTm[Ca2+]50refkUW[Ca2+]res[Ca2+]max τ CR τ CD
[][][][ μ M][ms 1 ][ μ M][ μ M][ms][ms]
Original3.22−1.215.470.5010.0460.0760.900164.4149.0
Mean3.55−1.214.350.5120.0710.0750.900142.4151.5
(33.6%)(0.5%)(23.6%)(2.3%)(59.3%)(0.1%)(2.8%)(28.3%)(10.4%)
SD1.21<0.011.130.0140.004<0.001<0.00152.318.8
(19.4%)(0.3%)(17.2%)(2.8%)(85.2%)(0.6%)(0.1%)(19.7%)(7.4%)
Min2.32−1.223.000.5000.0420.0750.89777.2123.3
(2.9%)(0.2%)(4.7%)(<0.1%)(1.1%)(0.2%)(<0.1%)(6.9%)(1.3%)
Max4.99−1.205.900.5380.1490.0770.900192.4178.8
(54.8%)(1.1%)(45.1%)(7.4%)(226.9%)(2.0%)(0.3%)(53.1%)(20.0%)
Table 10. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Pressure and volume biomarker values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Table 10. Results of the initial guess variation robustness analysis for the patient case 02-CoA. Pressure and volume biomarker values of the second step of the active mechanics personalisation approach are given. The original values produced by the default initial guesses are listed in the first row and means (M), standard deviations (SD), minima (Min), and maxima (Max) of five samples are listed in the subsequent rows. Means, standard deviations, minima, and maxima of the relative differences are given in brackets below.
Pmax dP dt | max dP dt | min PTD90SV
[mmHg][mmHg ms 1 ][mmHg ms 1 ][ms][mL]
Original1071.5−0.9446126
Mean1121.7−0.9422120
(6.6%)(12.8%)(5.5%)(9.7%)(10.7%)
SD90.30.14715
(7.4%)(15.6%)(3.3%)(6.6%)(6.5%)
Min1051.4−1.0351100
(0.4%)(0.2%)(0.8%)(3.6%)(1.5%)
Max1272.1−0.8473137
(19.4%)(41.6%)(9.5%)(21.3%)(20.6%)
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jung, A.; Gsell, M.A.F.; Augustin, C.M.; Plank, G. An Integrated Workflow for Building Digital Twins of Cardiac Electromechanics—A Multi-Fidelity Approach for Personalising Active Mechanics. Mathematics 2022, 10, 823. https://doi.org/10.3390/math10050823

AMA Style

Jung A, Gsell MAF, Augustin CM, Plank G. An Integrated Workflow for Building Digital Twins of Cardiac Electromechanics—A Multi-Fidelity Approach for Personalising Active Mechanics. Mathematics. 2022; 10(5):823. https://doi.org/10.3390/math10050823

Chicago/Turabian Style

Jung, Alexander, Matthias A. F. Gsell, Christoph M. Augustin, and Gernot Plank. 2022. "An Integrated Workflow for Building Digital Twins of Cardiac Electromechanics—A Multi-Fidelity Approach for Personalising Active Mechanics" Mathematics 10, no. 5: 823. https://doi.org/10.3390/math10050823

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop