Personalization of electro-mechanical models of the pressure-overloaded left ventricle: fitting of Windkessel-type afterload models

Computer models of left ventricular (LV) electro-mechanics (EM) show promise as a tool for assessing the impact of increased afterload upon LV performance. However, the identification of unique afterload model parameters and the personalization of EM LV models remains challenging due to significant clinical input uncertainties. Here, we personalized a virtual cohort of N = 17 EM LV models under pressure overload conditions. A global–local optimizer was developed to uniquely identify parameters of a three-element Windkessel (Wk3) afterload model. The sensitivity of Wk3 parameters to input uncertainty and of the EM LV model to Wk3 parameter uncertainty was analysed. The optimizer uniquely identified Wk3 parameters, and outputs of the personalized EM LV models showed close agreement with clinical data in all cases. Sensitivity analysis revealed a strong dependence of Wk3 parameters on input uncertainty. However, this had limited impact on outputs of EM LV models. A unique identification of Wk3 parameters from clinical data appears feasible, but it is sensitive to input uncertainty, thus depending on accurate invasive measurements. By contrast, the EM LV model outputs were less sensitive, with errors of less than 8.14% for input data errors of 10%, which is within the bounds of clinical data uncertainty. This article is part of the theme issue ‘Uncertainty quantification in cardiac and cardiovascular modelling and simulation’.


Methods (a) Data Pre-processing
In general, acquired clinical data are not suitable for being immediately interfaced with computational models, since data are not acquired synchronously at the same time. Rather, the time elapsed between acquisitions of the disparate data may be significant, in the range between minutes to days, and physiological conditions may be different, for instance, between acquisitions with and without sedation. These variabilities lead to inherent inconsistencies in the data which have to be considered in the process when identifying parameters of the underlying assumed model. Thus, careful pre-processing that factors in the specific circumstances of data acquisition is vital to obtain consistent input and validation data for modelling.
(i) Volumetric pre-processing Information on LV volume, V lv , is obtained from the clinical data in two distinct ways. The static LV volume at a given instant during a cardiac cycle was derived from 3DWH MRI scans acquired during diastasis. Segmentation of the static 3DWH image stacks yields the most accurate measure of the LV blood pool volume for one given instant in time. On the other hand, the change of LV volume over time and thus the flow in and out of the LV were obtained from SAX cine MRI scans which were segmented using the semi-automatic clinical analysis software syngo.via (Siemens Healthcare GmbH, Germany) following the guidelines proposed in [1]. Although carried out consecutively within a short time window <1 hour, discrepancies between the volumetric information obtained from these sources inevitably arise. SAX cine MRI scans were acquired under breath-hold conditions which increases intra-thoracic pressure that compresses the heart and, thus, entails a notable reduction in LV blood pool volume. On the other hand, 3DWH scans were acquired using a navigator under physiological breathing conditions. These differences in acquisition conditions, combined with additional confounding factors such as the use of different segmentation strategies for cine and 3DWH or changes in heart rate between acquisitions, lead to a systematic mismatch in the volumetric data. Since anatomical models of a given patient are built on the basis of 3DWH scans, cine MRI based volume traces must be adjusted to reconcile the volume of the anatomical model with dynamic volumetric data. In this study the 3DWH volume data were considered as a reference since these were acquired under more physiological conditions. Discrepancies between cine MRI and 3DWH MRI due to temporal misalignment as a consequence of altered heart rate were corrected by selecting a matching relative instant within a cardiac cycle based on where TT cine and TT 3DWH denote the MRI trigger times and T cine and T 3DWH the related cycle length of the analyzed heart beat. The difference between LV volumes obtained from the 3DWH segmentations (V 3DWH ) and from cine MRI corresponding to the cine trigger time (V TTcine ), is equivalent to the value of the horizontal shift ∆V of the volume curve (Eq. 1.2, see Fig. 1a). Thus, assuming that the dynamics of the cine MRI-based volume traces is accurate, but shifted by an offset relative to T 3DWH , volume traces can be reconciled by subtracting the offset For the sake of parameterizing the afterload model, the flow q across the aortic valve -which is q = − dV dt during ejection and q = 0 else (see Fig. 1a) -is the variable of primary interest, and not V lv per se. However, the rather poor temporal resolution of acquired volume traces -with SAX MRI approx. 30 ms to 40 ms are feasible -combined with additional jitter introduced by the segmentation procedures leads to oscillatory flow curves q(t) due to differentiation. These issues were addressed by smoothing the raw volume traces using cubic splines. Volume traces were split into two parts and interpolated from both ends towards each other. Polynomials of the third degree were used to enforce the derivative to meet dV dt = 0 at the onset of ejection, at end-systole and end-diastole.  (ii) Pressure pre-processing Pressure traces are recorded using a pig tail catheter measuring several pressure beats in the left ventricle p lv,m (solid blue line in Fig. 1b) as well as at the entrance port, arteria femoralis p lv,ref (dashed blue line), which is used as a reference. Consecutively, the catheter is pulled out of the ventricle through the aortic valve and pressure trace in aorta ascendens pao,m -shown in solid red together with the reference (dashed red) -is recorded. For the catheterisation procedure patient is under anaesthesia, whereas during echocardiography patient is in an awake state, which may lead to diverging measurement results (e.g. due to different HR). Thus, combining the outcome of both procedures for fitting should be avoided. Therefore, the pressure drop ∆pav,m across the aortic valve was calculated as the peak-to-peak difference of p lv,m and pm for CoA cases, whereas for AS cases results from echocardiography were used in combination with cuff measurements.
Recordings from catheterisation need to be pre-processed as during this highly invasive procedure systemic pressure may rise (visible as horizontal shift in reference traces) and HR may change. Further, even though measurements are ECG-triggered, individual traces do not start at the exact same instant of time within the cardiac cycle, making a time alignment necessary. Pressure traces contain several beats and the goal is to get a single representative averaged LV and aortic pressure beat (see Fig. 1b).
First the change in systemic pressure was corrected by finding the average minimum for all beats in both reference traces and by shifting pm and its reference p ref by the difference so that the two reference traces share the same baseline (see Fig. 1b).
To shift and scale the curves on the time axis some distinct points within the individual traces have to be identified, namely the starting point of the beats (aortic valve opening pressure pop,m), the peak valuespm (approximately minimum of p ref ), the closing pressure of the valve p cl,m and the end of the beat (equivalent to the start of the next beat). By simply looking for the starting and end point of a beat and scaling the beat as a whole, the resulting curve would not be representative for other heart rates. With decreasing heart rate only very small changes in duration of the ejection occur, whereas a significant difference is seen in the duration of diastole [2].
Further, the assumption was made that minima of pm beats are equivalent to the diastolic pressure and to pop,m. The time points of pop,m in pm corresponding to the value in p lv,m are used as starting points for the individual beats.
To find the dicrotic notch in pm, which equals the closing pressure of the aortic valve p cl,m , the time index of the peak reference pressure is used, followed by searching for a pressure minimum within a certain time frame (about 50 ms, window size may be varied depending on the shape of the dicrotic notch). The resulting p cl,m values need to be found in p lv,m to retrieve the corresponding time indices for valve closing. In this way, three time intervals (SYS1, SYS2 and DIA) which have to be aligned and scaled separately can be found to obtain pao,avg and p lv,avg (see Fig. 1b and Fig. 1c). The first two intervals make up the ejection phase (SYS) and the last interval corresponds to diastole (filling phase) of the cardiac cycle.
Additionally, RC-time was estimated with the decay time method [3] and solved via an optimization approach using the relation p(t) = α e − t τ for the decay of pressure during diastole (τ = RC). As from measured pao,m the time t ed and pressure p ed at the end of diastole is known where t 0 is the time point t cl of aortic valve closing plus 10% of the total cycle length Tcyc.

(iii) Pressure-volume synchronization
The averaged pressure beats pavg and p lv,avg were synchronized to match the cycle length of the volume traces. This was again done separately for the ejection phase (SYS) and the diastolic phase (DIA) visualized in Fig. 1c.

(b) Empirical Reference Data
To further corroborate physiological box constraints in AS cases for which no invasive data were recorded in the CARDIOPROOF study [4], a statistical analysis was conducted on a reference group of N = 290 patients treated for AS. All patients in the reference group had moderate to severe AS, undergoing invasive hemodynamic assessment to decide about the necessity of surgical or percutaneous treatment of the valve disease. During this procedure ventricular and aortic pressure are measured continuously with a 1.8 mm double-lumen fluid filled catheter in a resting condition in order to determine the transvalvular pressure gradient, ∆pav,m. The derived empirical values of the cohort (co) are summarized in Tab. 1 and are used as reference values for scaling. Table 1: Mean and standard deviation of empirical values obtained from evaluation of N = 290 pressure traces. Opening pressure pop,co of the aortic valve, closing pressure p cl,co and peak pressure in the aortapco are listed. MAP AUC,co is the mean arterial pressure estimated by the area under the curve (AUC) and RC l2,co is the diastolic decay time estimated by the decay time method in combination with an optimization approach using the l2-norm.

(c) Afterload Model
Values of ω i in Manuscript Eq. 2.3 are computed as follows: wherepco, pop,co and p cl,co result from the respective mean values of the cohort data (see Tab. 1) and the weighting factors γ i were chosen so that each term of the cost functional contributes to the same extent to the overall cost (γ 0 = γ1 = 100.0, γ 2 = 10.0).
The bounds b i used in Manuscript Eq. 2.4 were chosen to lie within physiologically realistic ranges and are derived from the cohort data (see Eq. 1.4). Their values are [5] was employed to generate electrical activation sequences which serve as a trigger for active stress generation in cardiac tissue. The hybrid R-E model combines a standard reaction-diffusion (R-D) model based on the monodomain equation with an eikonal model. Briefly, the eikonal equation is given as where (∇ X ) is the gradient with respect to the end-diastolic reference configuration Ω 0 ; ta is a positive function describing the wavefront arrival time at location X ∈ Ω 0 ; and t 0 are initial activations at locations Γ * 0 ⊆ Γ 0 with Γ 0 the boundary of Ω 0 . The symmetric positive definite 3 × 3 tensor V(X) holds the squared velocities (v f (X), vs(X), vn(X)) associated to the tissue eigenaxes, referred to as fibre (f 0 ), sheet (s 0 ), and sheet normal (n 0 ) orientations. The arrival time function ta(X) was subsequently used in a modified monodomain R-D model given as where an arrival time dependent foot current, I foot (ta), was added, which is designed to mimic subthreshold electrotonic currents to produce a physiological foot of the action potential. The key advantage of the R-E model is its ability to compute activation sequences at much coarser spatial resolutions that are not afflicted by the spatial undersampling artefacts leading to conduction slowing or even numerical conduction block, as it is observed in standard R-D models. Ventricular EP was represented by the Tusscher-Noble-Noble-Panfilov model of the human ventricular myocyte [6]. Note that activation sequences and electrical source distribution in the LV (1.5, 1.6) were computed in its end-diastolic configuration Ω 0 , that is, any effects of deformation upon electrotonic currents remained unaccounted for.
(ii) Active and passive mechanics in the LV and aorta The deformation of the heart is governed by imposed external loads such as pressure in the cavities or from surrounding tissue and active stresses intrinsically generated during contraction. Tissue properties of the LV myocardium and the aorta are characterized as a hyperelastic, nearly incompressible, anisotropic material with a non-linear stress-strain relationship. Mechanical deformation was described by Cauchy's equation of motion under stationary equilibrium assumptions leading to a quasi-static boundary value problem for t ∈ [0, T ], where u is the unknown displacement, F is the deformation gradient, S is the second Piola-Kirchhoff stress tensor, and (∇ X ·) denotes the divergence operator in the Lagrange reference configuration. The total stress S was additively decomposed according to where Spas and S act refer to the passive and active stresses, respectively. Passive stresses were modelled based on the constitutive equation given a hyper-elastic and transversely isotropic strain-energy function Ψ by Guccione et al. [7]. Here, the term in the exponent is and E = 1 2 (C − I) is the modified isochoric Green-Lagrange strain tensor, where C := J −2/3 C with J = det F. Default values of b f = 18.48, b t = 3.58, and b fs = 1.627 were used. The parameter C Guc was used to fit the LV model to an empirical Klotz relation [8] by a combined unloading and re-inflation procedure. In the aorta, unlike in previous studies [9], we refrained from assigning fibre structures, since our efforts were primarily focused on modelling the biomechanics of the LV and, to a lesser degree, the aorta.
Thus, in absence of information on structural anisotropy, an isotropic neo-Hookean model [10] was used The bulk modulus κ, which serves as a penalty parameter to enforce nearly incompressible material behaviour, was chosen as κ = 650 kPa in both (1.10) and (1.12).
A simplified phenomenological contractile model was used to represent active stress generation [11]. Owing to its small number of parameters and its direct relation to clinically measurable quantities such as peak pressure and the maximum rate of rise of pressure, this model is fairly easy to fit and thus very suitable for being used in clinical EM modelling studies. Briefly, the active stress transient is given by whereŜa is the peak isometric tension, φ(λ) is a nonlinear function dependent on fibre stretch λ = |Ff 0 | describing the length dependence of active stress generation, ts is the onset of contraction, τ C is the upstroke time constant, t dur is the active stress transient duration and τ R is the downstroke time constant. The active stress tensor in the reference configuration Ω 0 induced in fibre direction f 0 is defined as

(iii) Mechanics boundary conditions at the LV and aorta
The different BCs applied to the LV models are summarized in Fig. 2. The springs attached to the aortic rim and at the pericardium are shown in Fig. 2A, as well as the pressure BC in the cavity at the endocardial surface. The pericardial springs penalize displacement in normal direction only and are gradually scaled from the apex to the base. Therefore, the distance in apico-basal direction was used to create a penalty map, see Fig. 2A. To avoid non-physiological rotation, further springs were attached to the septum, see Fig. 2B. The location of the septal springs was selected automatically by constructing a local coordinate system spanned by the centers of the apical region, the Mitral valve (MV) and the Aortic valve (AV).

(e) Optimization
The optimization problem derived in Manuscript Sec. 2(b)i is non-convex and, hence, may possess several local minima. Thus, when solving this problem using an automated approach caution is warranted as the optimization may terminate with a non-competitive minimum. To circumvent this problem a combined global-local optimization method was applied. First, a neighborhood of a competitive local minimum was identified by stochastic sampling, which was then further refined by local optimization using adjoint-based derivatives. Specifically, stochastic Sobol sampling [12,13] was implemented to generate 5000 different parameter sets for Z, R and C within a sufficiently large box. The cost for each set was evaluated, the set of lowest cost determined and the box re-centered around this candidate set, see Fig. 3A. The procedure was iteratively repeated until the candidate set was found to lie within the center of the box, measured as 80% of the box size. The box centers are illustrated by dotted lines in Fig. 3A. Shrinking steps were applied then by incrementally narrowing down the search range for each parameter to better resolve the topography of the cost functional, see Fig. 3B.
The optimal result found from stochastic sampling was used then as an initial guess for local optimization. Here, quasi-Newton methods based on adjoint-based exact discrete derivatives were applied.

(a) Data Pre-processing
For N CoA = 7 cases, aortic pressure and volume traces were pre-processed according to Sec. 1(a)iiii in this supplement. Aortic valve opening pressure pop,m, peak pressurepm and valve closing pressure p cl,m were extracted from the processed pressure trace pm. End-diastolic LV pressure p lv,ed was obtained from the processed LV pressure trace p lv,m .
Due to the absence of measured pressure traces, only volume pre-processing, see Sec. 1(a)i, was applied to the N AS = 10 AS cases. Aortic valve opening pressure pop,m and peak pressure

(c) EM LV Model Parameterization (i) Anatomical modelling
Patient-specific anatomical models of all cases (N=17) were generated following established workflows mentioned previously and are visualized in Fig 4. (ii) Passive biomechanics Constitutive relations are represented in terms of the Guccione Model (see Eq. 1.10). From the anatomical models at ED state, a stress-free reference configuration is computed using default material parameters and p lv,ed in the LV (for AS cases the empirical value of 2.8 kPa is used, while for CoA cases patient-specific values extracted from invasive pressure measurements are taken). This is accomplished by unloading the geometry using a backward displacement method [14]. The EDPVR is fitted to the empiric Klotz EDPVR using default values for parameters b f = 18.48, b t = 3.58 and b fs = 1.627 from literature [7] adapting only the stiffness parameter C Guc . Values for the fitted passive material parameter C Guc are shown in Tab. 4 for all cases. Table 2: Pre-processed pre-treatment AS and CoA patient characteristics from SAX cine MRI, invasive catheterization, non-invasive cuff pressure recordings and echo-US including Sex, Age (y), heart rate (HR), end-diastolic volume (EDV), end-systolic volume (ESV), peak flow (q), catheter pressures (pop,m, p cl,m andpm) and the peak to peak pressure drop across the valve (∆pav,m) for CoA patient cases, diastolic and systolic cuff pressures (p dia,cuff and p sys,cuff ) and estimated closing pressure (p cl,m ) and pressure drop across aortic valve (∆pav,m) for AS patient cases. The impedance of the aortic valve Zv is computed as ∆pav,m overq.

(iii) Active stress model
The active stress model was fitted using default values of {τ C = 40 ms,Ŝa = 100 kPa, τ R = 110 ms} as initial guess. T dur was initialized with the RT interval observed in the ECG. A linear mapping was used to correct the active stress model parameters. In the AS cases where measured pressure traces p lv,m (t) were not available, onlyŜa was iteratively adjusted by a fixed-point iterationŜ a,i+1 =Ŝ a,i ·p lv /p lv,i . Rate of rise τ C and decline τ R in the active stress model were adjusted accordingly using the shape of the volume trace V lv,m . For the CoA cases, dp lv,m dt |max and dp lv,m dt | min were used additionally to fit τ C and τ R . As the focus of this study was on modelling systole, the diastolic isovolumetric relaxation and filling phases were not fitted. Fitted Parameters along with goodness of fit are summarized in Tab. 4 for all cases.

(d) Sensitivity Analysis
Sensitivity of fitting {Z, R, C} to uncertainty in input data {pop,m,pm, p cl,m } was analyzed for all N = 17 cases by independently varying input variables within a ±10% range. Relative sensitivities averaged pathology specific are summarized in Tab. 5. Further, the relation between relative deviation ot the Wk3 parameters resulting from variation in input data and the original input pressure ratio

Fitted Parameters
Goodness of Fit

Discussion (a) Limitations of the Study
For AS cases no invasive pressure recordings were available. Thus, systolic LV and central aortic pressures needed as input for the parametrization of the Wk3 model relied on a simplistic estimation where central aortic pressure was assumed to be equal to brachial systolic pressure   2.14 as measured by sphygmomanometry. It is well known that such estimates are highly inaccurate due to the large inter-individual variability in pressure wave augmentation [2,15]. While diastolic and mean arterial pressures are relatively constant throughout the arterial system and, thus, can be estimated from cuff measurements, this is not the case for systolic pressures. It is known that systolic pressure may be up to 40 mmHg (5.33 kPa) higher in the brachial artery than in the aorta [16]. This phenomenon of systolic pressure wave augmentation is attributed to an impedance mismatch between more compliant central vessels and the stiffer, less compliant more peripheral arterial vessels, which causes reflections of the forward travelling pressure wave. Central aortic pressure is therefore a composite of the forward travelling pressure wave and the reflected backward travelling wave, which combined augment pressure in the aorta relative to the brachial artery. How much the systolic pressure rises/diastolic pressure falls between aorta ascendens and the brachial artery is highly variable among patients depending on hemodynamic and pathophysiological characteristics such as body size and height, age, HR, stiffening of the aorta and general level of blood pressure [17]. This makes the extrapolation of pressure in the aorta ascendens from the sphygmomanometer measurements highly inaccurate. Beyond invasive measurements which may not always be viable non-invasive alternatives such as applanation tonometry exist, which may yield significantly more accurate estimates of systolic aortic blood pressure. For an overview of alternative non-invasive methods see, for instance, [18].
Beyond the significant measurement uncertainty, these wave transmission aspects are not taken into account by the Wk3 model. More elaborate Wk3 models that incorporate both reservoir-and wave-based aspects into a single lumped 0D model [19] or 1D models [20][21][22][23][24] may be able to better reflect vascular physiology, but have not been considered in this study. However, for providing an appropriate afterload to an EM model of the LV the use of a Wk3 model can be considered a fair trade-off. As shown in here, parameters can be identified fairly uniquely, the aortic pressure wave form approximates measurements sufficiently well and the physiological