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

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.


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 ([Ca 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, [Ca 2+ ] i , and fibre stretch, are generated and utilised to personalise the desired biophysically de-tailed 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.

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.

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].

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]  The spatiotemporal evolution of the displacement U in the LV and the AR is governed by Cauchy's first equation of motion that reads with initial conditions where S is the second Piola-Kirchhoff stress, ρ 0 is the density in the reference configuration, d 2 U dt 2 is the acceleration, and dU dt 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 Γ 0 R and Γ 0 N , respectively, with Γ 0 R ∪ Γ 0 N = ∂Ω 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 where is the strain energy function, whereC 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 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 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 Appendices A.1 and A.2. , 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 Robintype 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 = dV dt 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.

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 (Gradt a ) v(Gradt a ) = 1 in Ω 0 , (8) with initial conditions 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 with initial conditions 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 [Ca 2+ ] i evolution that provides [Ca 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 [Ca 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 [Ca 2+ ] i (electromechanical coupling) and the resulting deformation modifies the transmembrane potential (mechanoelectrical coupling). The latter can be caused through stretch-activated and stretch-modified ion currents and by stretch-induced changes of [Ca 2+ ] i . Stretch-induced changes of [Ca 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 Ca 2+ and from a modified sensitivity of the myofilaments to Ca 2+ . Electro-and mechano-mechanical coupling were considered but for the sake of simplicity, mechano-electrical coupling was not.

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 ms and 1.0 ms, 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 2022). 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. 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: where p LFM = {S max ref , τ SR ref , τ SD , t CR } are the parameters of the Tanh model (Appendix A.1) and B P = {P max , dP dt | max , dP dt | min , PTD 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: 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 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: 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, [Ca 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), [Ca 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.

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 [Ca 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 [Ca 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 ): with initial conditions λ(0) = λ 0 . (17) 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 ([Ca 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 [Ca 2+ ] i is used as input but [Ca 2+ ] i can also be generated from a coupled model of [Ca 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).
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: where t S amax = 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 j S a 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 [Ca 2+ ] i as input for the HFM, the HFM can also be coupled to a model of [Ca 2+ ] i evolution that is integrated in the chosen cellular EP model. This comes with the advantage that the model of [Ca 2+ ] i evolution can be included in the personalisation process. In this case, the constrained minimisation problem is extended to: where p CEM are the parameters of the [Ca 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 [Ca 2+ ] i trace. More specifically, CTD 50 and a CTD 90 -the [Ca 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 [Ca 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 , dS a dt | max , STD 30 , STD 50 , STD 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 Tables A2-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.
Algorithm 1 Two-step multi-fidelity approach for personalising biophysically detailed active mechanics models. Compute error between simulated (B sim p ) and clinically measured (B clin p ) pressure biomarkers.

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 S1) or by the fraction of variance attributed to each model input including interactions with the other model inputs (total-effect sensitivity index ST). 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).
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 organscale 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].

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.

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. ] 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.

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 , dP dt | max , dP dt | min , PTD 90 , SV). 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 PTD 90 had means (standard deviations; SD) of 5.4% (SD: 3.4%) and 1.4% (SD: 2.2%), respectively. For dP dt | max and dP dt | 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 SV 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   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 PTD 90 remained almost unchanged: 5.5% (SD: 2.3%) and 2.9% (SD: 1.5%), respectively. For dP dt | max and dP dt | 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 SV also remained almost unchanged 11.1% (SD: 6.2%).

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. Tables 3-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 dP dt | max (10.5%; SD: 7.3%). The pressure-volume loops and the individual pressure and volume traces are compared in Figure 6.

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. Tables 7-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 max ref , τ SD , T CR , and the VFM parameter R MVf were below 2.9%. For τ SR 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  ] 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. 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 (Sim 1 , Sim 2 ) 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. 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.

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 [Ca 2+ ] i evolution on a purely phenomenological basis. This was done because information on [Ca 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].

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 [Ca 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).

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 τ SR 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 , [Ca 2+ ] 50 ref ), two Rice model parameters ([Ca 2+ ] res , [Ca 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].

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.

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 noninvasive 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 [Ca 2+ ] i evolution was simplified. There are numerous biophysically detailed models of cellular EP and [Ca 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.

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. 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: The Tanh model [37,38] describes the evolution of cellular active stress as function of fibre stretch and the electrical activation time: Here, S max ref is the maximum isometric cellular active stress and 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, 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. The Land model [39] describes the evolutionn of cellular active stress as function of [Ca 2+ ] i and both the fibre stretch λ and the fibre stretch rate dλ dt . 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 Ca 2+ , troponin C, troponin I, and tropomyosin that control the availability of myosin binding sites on actin. The dynamics of the interaction between Ca 2+ and troponin C is described by where CaTRPN is the fraction of regulatory troponin C sites with bound Ca 2+ , k TRPN represents the unbinding rate of Ca 2+ from troponin C, and n TRPN is the cooperativity of the binding between Ca 2+ and troponin C.
where B is the fraction of blocked myosin binding sites on actin, and k U are the troponin I and tropomyosin rate constants, respectively, n Tm is the steadystate relation between CaTRPN 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 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 with the steady-state ratios The distortion-depending unbinding rates of the cross bridges are and these are coupled to a distortion-decay model given by Here, ξ W and ξ S are the stretch rate-dependent mean distortions and are related to the magnitude of the instantaneous response to the distortion with some scaling A eff , whereas 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 where S max ref is the maximum isometric active stress and fibre stretch effects on the generated active stress are phenomenologically captured by where β 0 represents the change in maximum active stress based on changes in filament overlap. The parameter values are given in Table A2.

. Land Model of Cellular Passive Stress Evolution
The Land model [39] describes the evolution of cellular passive stress by a threeelement 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 with the stress components the series strain constraint and the series stress constraint 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.
with  [43]. To this end, the unconstrained minimisation problem 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.

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 [Ca 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.

Appendix C. Biomarker Definitions
The definitions of the biomarkers used in this study are illustrated in Figure A1. Based on the cavity pressure trace, the pressure biomarkers are defined as: Here, P max is the maximum, dP dt | max and dP dt | min are the maximum and minimum rates, and PTD 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: Here, SV 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) SV (A27) against t from t 4 to t 5 and the slope s VR is computed by an ordinary linear regression of against t from t 6 to t 7 .
Based on the cellular active stress trace, the active stress biomarkers are defined as: dS a dt , STD 30 = t 8 − t 12 , STD 50 = t 9 − t 13 , STD 90 = t 10 − t 14 . (A29) Here, S a max is the maximum, dS a dt | max is the maximum rate, and STD 30 , STD 50 and STD 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: (A31) 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.  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.