## Abstract

We introduce a parameter estimation framework for automatically and robustly personalizing aortic haemodynamic computations from four-dimensional magnetic resonance imaging data. The framework is based on a reduced-order multiscale fluid–structure interaction blood flow model, and on two calibration procedures. First, Windkessel parameters of the outlet boundary conditions are personalized by solving a system of nonlinear equations. Second, the regional mechanical wall properties of the aorta are personalized by employing a nonlinear least-squares minimization method. The two calibration procedures are run sequentially and iteratively until both procedures have converged. The parameter estimation framework was successfully evaluated on 15 datasets from patients with aortic valve disease. On average, only 1.27 ± 0.96 and 7.07 ± 1.44 iterations were required to personalize the outlet boundary conditions and the regional mechanical wall properties, respectively. Overall, the computational model was in close agreement with the clinical measurements used as objectives (pressures, flow rates, cross-sectional areas), with a maximum error of less than 1%. Given its level of automation, robustness and the short execution time (6.2 ± 1.2 min on a standard hardware configuration), the framework is potentially well suited for a clinical setting.

## 1. Background

Blood flow computations, when used in conjunction with patient-specific anatomical models extracted from medical images, provide important insights into the structure and function of the cardiovascular system. These techniques have been proposed for diagnosis, risk stratification and surgical planning [1]. For accurate patient-specific computations, the role of physiologically sound boundary conditions (inlet and outlet boundary conditions, and vascular wall properties) is well appreciated in the literature.

Various approaches have been proposed for personalizing the inlet and outlet boundary conditions. The outlet boundary conditions model the effect of the distal vasculature and are typically represented by lumped parameter models. The most widely used lumped parameter model is the three-element Windkessel model [2].

Various calibration procedures for the outlet boundary conditions have been proposed based on multiple pressure and/or flow measurements, which are typically available in a clinical scenario. A fully automatic optimization-based calibration method for the Windkessel models was suggested [3]. This method was then further refined to both increase its robustness as well as to reduce the number of iterations required for reaching convergence [4]. Furthermore, an adjoint-based method for calibrating the Windkessel parameters was proposed [5]. A competitive alternative to the above mentioned optimization-based methods is represented by filtering-based methods [6].

The vascular wall properties determine the arterial distensibility, which is an important factor for the development and assessment of cardiovascular diseases [7]. Typically, the arterial distensibility is described by the arterial compliance (or the arterial elastance—the inverse of the compliance), which is responsible for important functional aspects of the systemic circulation: larger blood flow rate in the coronary arteries during diastole, reduction of left ventricular afterload (during systole) and continuous flow at the level of the capillaries. Previous studies indicate that arterial compliance changes with age [8] and hypertension [9].

The arterial wall properties at a certain location in the systemic circulation can be described by the local compliance, typically defined as area compliance, *C*_{A}, or by the local pulse wave velocity *c*. The pulse wave velocity (PWV) is also used as a robust prognostic parameter in preventive cardiovascular therapy [10]. Alternatively, the arterial wall properties can be described globally, for a certain region or for the entire systemic circulation, by the volumetric compliance *C*_{V}.

Several approaches have been proposed in the past for non-invasive estimation of arterial wall properties [11]. Many of them rely on the transit time of the flow/pressure wave, i.e. the time that a flow/pressure wave needs to travel the distance between two locations. These methods may have a low accuracy if the distance used for the estimation of the transit time is relatively short, and can only provide an average value of the regional mechanical wall properties for the region of interest. Other methods estimate local vascular wall properties:

— the ACM method [12] estimates the local area compliance as

*C*_{A}= Δ*A/*Δ*P*, where Δ*A*is the difference between minimum and maximum cross-sectional area during a heart cycle, and Δ*P*is the pulse pressure. The wall properties are estimated independently at each location and do not take into account the overall haemodynamic flow conditions;— the pulse pressure method (PPM) [13] estimates the downstream volumetric compliance from the flow rate waveform and from the pulse pressure. A single value is estimated for an entire vascular segment;

— the PU-loop method [14] estimates the local PWV as

*c**=*d*P*/*ρ*d*U*during early systole, where d*P*is the derivative of the pressure and d*U*is the derivative of the blood flow velocity. The use of derivatives in the formulation amplifies the measurement/computation errors;— the DU-loop method [15] estimates the local PWV as

*c**=*0.5 d*U*/d(ln*D*) during early systole, where d*U*is defined as above and*D*is the diameter. Similar to the PU-loop method, the use of derivatives in the formulation amplifies the measurement/computation errors.

Furthermore, previously reported patient-specific blood flow computation frameworks typically set the vascular wall properties either based on a best fit to experimental data [16], or determine a single PWV value for the entire anatomical model [17,18].

In this paper, we propose a parameter estimation framework for automatically and robustly personalizing aortic haemodynamic computations from four-dimensional magnetic resonance imaging (4D MRI) data. The framework is based on a reduced-order multiscale fluid–structure interaction (FSI) blood flow model and personalization procedures. The latter calibrate inlet and outlet boundary conditions, as well as the regional mechanical wall properties, to ensure that the computational results match the patient-specific measurements.

## 2. Material and methods

### 2.1. Extraction of four-dimensional anatomical and flow information from medical imaging data

The parameter estimation and computational framework introduced herein is based on anatomical and flow information extracted from 4D flow MRI medical imaging data. A proprietery tool (research prototype, not for diagnostic use) is used to reconstruct and pre-processes the data by applying a series of image correction algorithms such as phase anti-aliasing and motion tracking (requiring 2–3 min of processing time) [19].

Once the 4D flow data are prepared, the steps displayed in figure 1 are performed sequentially. First, the user initiates a semi-automatic segmentation procedure by selecting a set of seed points, starting at the inlet of the ascending aorta and ending at the outlet of the descending aorta. Typically two to four seed points are required for a reliable initialization, and a preliminary, rough segmentation is performed based on a clustering approach that groups voxels in the 4D flow data into static tissue, air/lung and blood. Next, the centreline of the aorta is automatically extracted from it. Once the user has verified the proposed centreline, the segmentation is automatically refined, and the refined geometry may be manually corrected if deemed appropriate. As a result, the vessel lumen anatomical model and its centreline are obtained. Subsequently, a large number (typically 50) of cross-sectional planes (analysis planes) along the aorta are automatically generated. For each analysis plane, the flow data are systematically analysed to determine the time-varying cross-sectional areas and the time-varying flow rates along the aorta. Anatomical and flow information is thus extracted only for the aorta (ascending aorta, aortic arch and descending aorta), and not for the supra-aortic branches.

In addition to flow and geometrical data, we further rely on non-invasive, cuff-based pressure measurements obtained from the left arm.

### 2.2. Reduced-order multiscale fluid–structure interaction blood flow model

The parameter estimation framework for personalizing aortic haemodynamic computations employs a reduced-order multiscale FSI blood flow model, which is based on a quasi-one-dimensional and a lumped parameter model (zero-dimensional model) [17].

The one-dimensional blood flow model is derived from the three-dimensional Navier–Stokes equations based on a series of simplifying assumptions. The resulting governing equations ensure mass and momentum conservation [20]. A state equation, which relates the pressure inside the vessel to the cross-sectional area, is used to close the system of equations. The vessel wall is modelled as an elastic material:
2.1where *x* denotes the axial location, *t* denotes the time, *A*(*x*, *t*) is the cross-sectional area, *p*(*x*, *t*) is the pressure, *E* is the Young modulus, *h* is the wall thickness, *r*_{0} is the initial radius corresponding to the initial pressure *p*_{0} and *A*_{0} is the initial cross-sectional area. At each bifurcation, the continuity of flow and total pressure is imposed.

In the past, we have also employed viscoelastic wall laws, indicating that viscoelasticity leads to changes of up to 6% in the pressure values and up to 20% in the flow rate values in a full body arterial model [20]. Another effect of wall viscoelasticity is the more pronounced dampening of high-frequency oscillations observed in the computed time-varying pressures and flow rates [21]. However, these observations apply predominantly to the distal locations of the human arterial system, which are not modelled herein.

Time-varying flow rate profiles are used as inlet boundary condition, while three-element Windkessel models are coupled at the outlets of the one-dimensional model.

Such reduced-order models have been shown to accurately predict time-varying flow rate and pressure wave forms [22]. Recent research activities have shown the growing interest in the reduced-order blood flow models for specific parts of the circulation under patient-specific and pathologic conditions: the coronary circulation [23], the abdominal aorta [24], proximal part of the aorta [17], the aortic valve [25] and the peripheral circulation [26]. The reduced-order model used in this study has been previously introduced in [17,23].

### 2.3. Parameter estimation framework

Figure 2 displays an overview of the proposed parameter estimation framework for personalizing the aortic haemodynamic computations. The reduced-order multiscale FSI blood flow model is first initialized and two independent calibration procedures are sequentially and iteratively employed for automatically and robustly personalizing the aortic haemodynamic computations:

Personalization of the Windkessel parameters in the outlet boundary conditions of the multiscale circulation model.

Personalization of the regional mechanical wall properties of the aorta.

The two calibration procedures are run until the convergence criteria are met for both of them simultaneously. Finally, the results of the haemodynamic computations are post-processed to determine clinically relevant characteristics.

#### 2.3.1. Initialization of the reduced-order multiscale fluid–structure interaction blood flow model

Figure 3*a* displays the multiscale FSI blood flow model for the aorta and the supra-aortic branches. The aorta is divided into multiple segments for which the regional mechanical properties are estimated (these are numbered from 1 to 7 in figure 3*a*, but the actual number of segments can vary and depends on the length of the ascending and descending aorta).

The initialization of the blood flow model consists of the following five steps, which are described in detail below:

(1) Defining the bifurcation locations of the supra-aortic branches.

(2) Defining the average flow rate values for the ascending and descending aorta.

(3) Defining the one-dimensional segments and their geometry.

(4) Defining the inlet boundary condition and the initial parameter values at the outlet boundary condition.

(5) Defining the initial regional mechanical wall properties.

As was mentioned in §2.1, anatomical and flow rate information is only available for the aorta. However, because the supra-aortic branches draw away from the aorta a significant volume of blood (30–50%), to run accurate haemodynamic computations we augment the anatomical model with supra-aortic branches. To this end, the bifurcation point of the first supra-aortic branch (the brachiocephalic artery) is determined by navigating downstream, starting from the ascending aorta inlet, through the centreline locations until a cross-sectional plane is found, for which the average flow rate decreases below 90% of the average flow rate at the upstream centreline locations, and no downstream cross-sectional plane with a larger flow rate exists. Similarly, the bifurcation point of the third supra-aortic branch (the left subclavian artery) is determined by navigating upstream, starting from the descending aorta outlet, through the centreline locations until a cross-sectional plane is found for which the average flow rate is larger than 110% of the average flow rate at the downstream centreline locations, and no upstream cross-sectional plane with a lower flow rate exists. Finally, the bifurcation point of the second supra-aortic branch (the left common carotid artery) is set midway between the other two supra-aortic branches.

Next, average flow rate values are estimated for the ascending and the descending aorta. Owing to measurement noise, the average flow rate varies slightly between consecutive cross-sectional locations. A linear least-squares fit-based algorithm is employed, which is used to filter out locations with very large or very low average flow rate values. Based on the remaining locations, a final average flow rate value is determined for the ascending and the descending aorta.

In the following, the one-dimensional segments and their geometry are defined. First the number of segments is set for each branch: two segments for the aortic arch, and multiple segments for the ascending (at least two) and descending aorta (at least three). Spatially varying cross-sectional area values are defined for each segment to obtain a geometry which is reliably reflecting the actual three-dimensional geometry. The initial cross-sectional area values are based on the end-diastolic phase.

Population-average geometric properties [22], which are scaled based on the patient-specific aorta size, are applied to define the one-dimensional segments corresponding to the supra-aortic branches (a fixed length of 2 cm is used for each supra-aortic branch).

Subsequently, the inlet boundary condition is defined: the flow rate profile measured at the first analysis plane is scaled so as to match the average ascending aorta flow rate value estimated as described above. Next, three initial parameter values need to be specified at each outlet. First we compute the average pressure at the inlet of the left subclavian artery, following an approach validated in [12]. We use as input the brachial systolic pressure *P*_{b-s} and the brachial diastolic pressure *P*_{b-d}. The diastolic pressure at the inlet of the left subclavian artery is set equal to *P*_{b-d} (pulse wave propagation phenomena affect the systolic pressure but not the end-diastolic pressure, because at end diastole the flow in the large arteries is approximately zero, and wave propagation effects are negligible), while the systolic pressure is computed from
2.2

Next, the mean arterial pressure at the inlet of the left subclavian artery is computed using a form factor of 0.4: 2.3

Since the variation of the average arterial pressure in the aorta is typically small, is used for determining the initial total resistance at each outlet as ratio between average pressure and average flow rate: 2.4

To determine the average flow rate for each supra-aortic branch the total supra-aortic average flow rate, *Q*_{supra-aortic}, is computed as difference between the average flow rates in the ascending and the descending aorta. This flow is then distributed to the branching vessels proportionally to the square of the radius. To minimize reflections, the proximal resistance of each Windkessel model is set equal to the characteristic resistance of the corresponding outlet segment, while the distal resistance is computed as difference between total and proximal resistance.

For the estimation of compliance values, we first consider a population-average compliance value (*C*_{tot}) [24], which is then distributed to the four outlets proportionally to the square of the radius.

To initialize the regional mechanical wall properties we first determine the arterial wall properties at the bifurcation of the left subclavian artery. Equation (2.1) is rewritten, based on *P*_{LSA-s} and *P*_{LSA-d}, as
2.5where *A*_{LSA-s} and *A*_{LSA-d} are the maximum (systolic) and minimum (diastolic) cross-sectional area values, and *β* represents the wall stiffness. Hence
2.6

This value is used to initialize the stiffness for the entire aorta. To estimate the wall properties of the supra-aortic vessels, we use a slightly modified approach, under which the wall properties of each supra-aortic segment are computed separately. This is done to minimize the wave reflections at the bifurcations [17]. The initial pressure in the entire domain is set equal to *P*_{b-d}.

#### 2.3.2. Parameter estimation procedure for personalizing the outlet boundary conditions

The objective of the parameter estimation procedure described in this section is to adapt the parameters of the Windkessel models coupled to the outlets of the one-dimensional model, under the constraint that the blood flow solutions should (i) maintain the same average flow-split at each outlet as determined with the procedure described in §2.3.1 and (ii) replicate the measured systolic and diastolic pressure at the inlet of the left subclavian artery. Out of the four flow-split values (three supra-aortic branches and the descending aorta) only three are used as objectives, because the fourth one is obtained automatically as difference.

The parameter estimation problem is formulated as a solution to a system of nonlinear equations, with each equation representing the residual error between the computed and measured quantity of interest. To determine the values of all the residuals (** f**(

*x**)), a computation with the parameter values*

_{i}

*x**is required. Since the absolute values of the adapted parameters and of the residuals generally differ by orders of magnitude, for the calibration method both the parameter and the objective residuals have been scaled using typical values.*

_{i}The parameters to be estimated are the total resistances of the three supra-aortic vessels and of the descending aorta, and the total compliance. The following system of nonlinear equations is numerically solved to obtain the optimum value of each parameter:
2.7where *P*_{max-LSA} is the maximum (systolic) pressure, *P*_{min-LSA} is the minimum (diastolic) pressure, both at the inlet of the left subclavian artery, represents a flow rate split, while (•)_{comp} refers to a value computed using the lumped parameter/multiscale model and (•)_{ref} refers to the reference value. Index LCC refers to the left common carotid artery, LS to the left subclavian artery and DAo to the descending aorta.

The typical values of the objectives are set equal to 1 mmHg for the pressure-based objectives and to *Φ*/100 for the flow rate split-based objectives.

The nonlinear system of equations is first solved for a lumped parameter model, composed of the Windkessel models used in the multiscale model (figure 3*b*). The initial solution *x*_{0} is determined using the steps described in the previous section and a dogleg trust region algorithm is applied to iteratively determine the solution of the nonlinear system of equations. The solution determined for the lumped parameter model is then adapted as described in [5], to compensate for the haemodynamic properties (resistance and compliance) of the multiscale model that are not taken into account in the lumped parameter model. As a result, the risk of a failure in finding a solution to the nonlinear system of equations is reduced, and, importantly, the number of calibration iterations required to reach the final solution is decreased.

Next, equation (2.7) is solved for the multiscale reduced-order blood flow model. Each computation, with a given set of parameter values, is run until the *L*_{2} norms of the normalized differences between the pressure and flow rate profiles at the current and the previous cardiac cycle are smaller than 10^{−5}. A quasi-Newton method is employed at this stage, where at each iteration the Jacobian is only updated and not recomputed, to ensure short computation times. If all objective residuals are smaller than the tolerance limit (taken here equal to ), the calibration method is terminated.

At each iteration, the parameter values of the three-element Windkessel models (figure 3) are set as follows, based on the parameter values in equation (2.7): the proximal resistance of each Windkessel model is set equal to the characteristic resistance of the corresponding outlet segment, while the distal resistance is computed as difference between total and proximal resistance. The total compliance is distributed to the four outlets proportionally to the squares of their respective radii.

#### 2.3.3. Parameter estimation procedure for personalizing the regional mechanical wall properties

The objective of the parameter estimation procedure described in this section is to adapt the local wall stiffness along the aorta so as to obtain a good match between the measured and the computed cross-sectional area variation at the analysis planes. The parameter estimation procedure is based on a nonlinear least-squares method, which minimizes the following cost function:
2.8where *m* is the total number of measurements, *j* refers to a specific measurement location along the aorta and *r _{j}*(

**x**) are the residuals computed as difference between the measured and the computed quantities: 2.9where is the measured maximum variation in the cross-sectional area during a heart cycle at location

*j*, and is the computed maximum variation in the cross-sectional area during a heart cycle at location

*j*.

The parameter vector **x** contains the wall stiffness at the start and end of each one-dimensional segment in the computational model (a linearly varying stiffness is imposed between the inlet and outlet of each one-dimensional segment):
2.10

The cost function is minimized using a quasi-Gauss–Newton method (a modified Newton's method with line search). Thus, we first determine a search direction , and choose a step length *α _{k}* to satisfy the Armijo and Wolfe conditions [27]. At each iteration, once the wall properties of the aorta have been updated, the wall properties of the supra-aortic branches are modified to minimize the wave reflections at the bifurcations (see §2.3.1). Similar to the set-up described in the previous section, the Jacobian is only computed once, and then updated at each further iteration. Once the cost function converges (its variation from one iteration to the next becomes smaller than 1%), the calibration procedure is terminated.

As depicted in figure 2, once the calibration procedure for determining the regional mechanical wall properties has converged, the convergence criteria of the calibration procedure for determining the outlet boundary conditions are verified. If these are not satisfied, the two parameter estimation procedures are rerun. The convergence criteria of the first calibration procedure may no longer be satisfied once the second parameter estimation procedure has been applied, because a change in the wall properties generally induces a change in the pressure and flow rate values.

#### 2.3.4. Computation of clinically relevant characteristics

Once the calibration procedures of the parameter estimation framework have converged, the computational results are post-processed to determine a series of clinically relevant characteristics.

The local pulse wave velocity: 2.11

The local area compliance:

a. Using the wall stiffness: 2.12

where

*A*(x) is the average value of the cross-sectional area at location*x*.b. Using the ACM method: 2.13

where Δ

*A*(*x*) is the maximum variation of the cross-sectional area at location*x*and PP(x) is the pulse pressure at location*x*, determined from the haemodynamic computations.

The downstream volumetric compliance: 2.14

where PPM refers to the pulse pressure method [28], which uses as input the time-varying flow rate at location

*x*and the pulse pressure at location*x*, as determined from the haemodynamic computations.

## 3. Results

To evaluate the performance of the proposed parameter estimation framework, we present results for 15 datasets acquired from patients with aortic valve disease (60% male) enrolled in the ‘Cardioproof’ project [29]. The patients' ages ranged from 9 to 69 years (mean: 30.9 ± 21.8 years). The study complied with the Declaration of Helsinki for investigation in human beings. The study protocol was approved by the local ethics committee and each patient signed an informed consent form before the enrolment in the study.

Blood was modelled as an incompressible Newtonian fluid with a density of 1.050 g cm^{−3} and a dynamic viscosity of 0.040 dyn cm^{−2} s^{−1}. The grid size for the numerical solution of the one-dimensional model was 0.05 cm, while the time-step (limited by the CFL-condition) was set equal to 2.5 × 10^{−5} s.

### 3.1. Personalization of the haemodynamics computations

The calibration procedures converged successfully for all 15 datasets, with an average execution time of 6.2 ± 1.2 min on an off-the-shelf computer (Intel i7 processor). Table 1 displays the calibration results: the objective and parameter values of the outlet boundary conditions calibration procedure (reference objective values are displayed in parenthesis), and, for the regional mechanical wall properties calibration procedure, the initial and final values of the cost function (equation (2.8)) are displayed alongside the estimated wall stiffness. The number of calibration iterations is displayed both at global level (iterations of the loop in figure 2), and separately for each calibration procedure (as a sum of all calibration iterations during the different global iterations). Only one or two global iterations were required to reach convergence.

The five objectives defined in (2.7) are matched closely for each dataset and the total number of calibration iterations for personalizing the outlet boundary conditions was at most three (1.27 ± 0.96). A very good match between the reference and the computed values of the objectives was observed (mean relative difference < 0.5%).

The total number of calibration iterations for personalizing the regional mechanical wall properties was 7.07 ± 1.44; the value of the cost function was decreased from 2.95 ± 1.90 to 0.23 ± 0.37. The mean wall stiffness for the entire cohort was (0.93 ± 0.37) × 10^{3} mmHg. This value is similar to the mean wall stiffness of large arteries described previously in the literature, obtained by performing a best fit to experimental data (0.85 × 10^{3} mmHg) [16].

Next, we present in figure 4 as an example the calibration results obtained by personalizing the boundary conditions of patient dataset 15: both objective and parameter values, obtained by running the parameter estimation procedure described in §2.3.2, are displayed (the dotted red lines represent the reference values of the objectives). Two global iterations were required to reach final convergence, and during each global iteration, one local iteration was required to match the objectives related to the outlet boundary conditions (2.7).

Iteration 0 refers to the results obtained by running the computation with the parameter values determined by solving the nonlinear system of equations for the lumped parameter model, and by applying the algorithms for compensating for the haemodynamic properties (resistance and compliance) of the multiscale model. These initial parameter values lead to a good initial match between target and actual values of the objective. In particular, the computed values of *Φ*_{LCC}, *Φ*_{LS} and *Φ*_{DAo} are within 1% of their target values. *Φ*_{LCC}, *Φ*_{LS} and *Φ*_{DAo} are less affected when the algorithm switches from the lumped parameter model to the multiscale model than the maximum and minimum pressures. This is due to the fact that the former are mainly determined by the resistances, whereas the maximum and minimum pressures are influenced by both resistances and compliances. In turn, the total resistance of the model is mainly determined by the distal vasculature, whereas the compliance is significantly influenced (increased) by the proximal aorta model.

The objective values at the end of the first global iterations are different from the objective values at the start of the second global iteration. This is due to the fact that the regional mechanical wall properties are calibrated in between, thus also affecting the objectives related to the outlet boundary conditions, and leading to the requirement of running a second global iteration. The parameter values, on the other hand, are the same at the end of the first global iteration and at the start of the second global iteration.

Next, we present in figure 5 the calibration results obtained by personalizing the regional mechanical wall properties of patient dataset 15: both objective and parameter values, obtained by running the parameter estimation procedure described in §2.3.3, are displayed. Four and two local iterations were required during the first and second global iteration, respectively.

The cost function value at the end of the first global iterations is slightly different from the value at the start of the second global iteration, because the parameters of the outlet boundary conditions are calibrated in between. Initially, the wall stiffness is identical at all locations (as determined during the initialization of the multiscale model (2.6)), but, as the calibration procedure progresses, the wall stiffness of each segment is personalized so as to match the measured cross-sectional area variation.

In figure 6*a*, we display a comparison of measured (4D MRI) and computed minimum and maximum cross-sectional areas along the centreline of the aorta for patient dataset 15: the values match closely, indicating that the parameter estimation procedure calibrating the wall properties is able to provide a reliable personalization of the arterial wall stiffness. The small differences that can be observed are partially due to the approximation of the three-dimensional geometry through the one-dimensional model. Furthermore, figure 6*b* and figure 6*c* display a comparison of the computed and measured time-varying cross-sectional areas at one ascending and one descending aorta location. Importantly, not only the absolute values, but also the timing of the time-varying profiles match closely, confirming that the wave propagation phenomena are captured correctly by the personalized haemodynamic model.

### 3.2. Regional mechanical wall properties

Table 2 displays the quantities of interest related to the regional mechanical wall properties of the 15 patient datasets included in the study, computed as described in §2.3.4. Alongside the PWV determined from (2.11) we also included PWV estimates determined with the transit time method [30]: *c* = Δ*x/*Δ*t*, where Δ*x* is the length along the centreline of the aorta and Δ*t* is the time required for the flow rate wave to travel along this distance. Since Δ*t* is typically very small we have divided the aorta into only two regions and the transit time based PWV was estimated separately for these two regions. The first region contains the ascending aorta and the aortic arch, while the second region contains the descending aorta. The time Δ*t* is computed as the interval between the onset (foot) of the flow curves at the start and end of a region. The location of the onset (foot) is determined by the intersection point of the upslope curve and the minimum flow rate. The upslope curve is approximated by the line connecting the points at 30% and 70% of the maximum flow rate at the particular location.

As can be observed in table 2 the mean PWV estimated using (2.11) lies typically in between the two values estimated with the transit time method. Owing to the small values of Δ*t*, the transit time method is very sensitive to measurement errors and noise. As an example, we display in figure 7 a comparison of the local PWV and the transit time based PWV for patient dataset 15. The local PWV is indeed slightly larger in region 1 than in region 2, but the transit time method overestimates its value in the first region, and underestimates it in the second region.

Furthermore, table 2 displays the local area compliance, as determined from the wall stiffness and using the ACM method. The two methods lead to similar area compliance results, with a mean difference of only −0.359 × 10^{−3} ± 0.166 × 10^{−3} cm^{2} mmHg^{−1}.

Finally, table 2 also displays the downstream volumetric compliances as determined at the inlet of the ascending aorta and the inlet of the descending aorta. The average total volumetric compliance is 988.4 × 10^{−6} ± 297.7 × 10^{−6} cm^{4} s^{2} g^{−1}, which is in the range of previously reported values in the literature. On average, the descending aorta volumetric compliance is 29% smaller than the ascending aorta volumetric compliance.

Figure 8 displays for patient dataset 15 a comparison of the local area compliance, as determined from the wall stiffness and using the ACM method, and the downstream volumetric compliance along the centreline of the aorta. The regions with large local area compliance are the regions with low local PWV (figure 7) and vice versa (the area compliance and the PWV are inversely proportional). The downstream volumetric compliance decreases significantly towards the descending aorta. A marked decrease can be observed once the bifurcation points of the supra-aortic branches have been passed. This is due to the fact that the downstream compliance at an ascending aorta location contains the volumetric compliance of all arteries supplying the arms and the cerebral circulation.

## 4. Discussion and conclusion

This paper introduced a parameter estimation framework for automatically and robustly personalizing aortic haemodynamic computations from 4D MRI medical image data. The haemodynamic computations are based on a reduced-order multiscale FSI blood flow model, which is fully personalized: the patient-specific anatomical model and the time-varying flow rate inlet boundary condition are derived directly from the medical imaging data, and the proposed parameter estimation framework calibrates the parameters of the outlet boundary conditions and the regional mechanical wall properties along the centreline of the aorta so as to match the clinical measurements. The former are personalized by solving a nonlinear system of equations, while the latter are personalized by minimizing a cost function. Owing to the different nature of the calibration procedures, these are run sequentially and iteratively until both of them are converged. To the best of our knowledge, this is the first time that spatially varying regional mechanical wall properties are personalized in the context of patient-specific haemodynamic computations.

We evaluated the methodology by investigating 15 patient datasets. Since all of them converged successfully within only two global iterations, three iterations for the calibration of the outlet boundary conditions and 10 iterations for the calibration of the regional mechanical wall properties, the proposed method is deemed robust. The estimated aortic pressures and flow rate distributions between supra-aortic branches and descending aorta were matched with an error of less than 1%. We used three-element Windkessel models as outlet boundary conditions, but other types of physiologically sound boundary conditions, like the structured tree boundary condition, may be employed instead [31]. The cost function related to the estimation of regional mechanical wall properties, computed based on the differences between the cross-sectional area variations in the 4D MRI data and the blood flow model, was minimized to an average value of 0.23 cm^{2} for the entire cohort. Moreover, as can be observed in figure 6, not only the cross-sectional area variation (difference between minimum and maximum value) is matched well, but also the time-varying cross-sectional area profiles typically match the measured profiles.

In addition to its robustness and excellent capability of matching the clinical measurements, the proposed method is computationally also very efficient: the average execution time was only 6.2 ± 1.2 min on a standard hardware configuration. There are several aspects contributing to this achievement. First of all, we use a computationally efficient reduced-order FSI model. Secondly, the parameters in the outlet boundary conditions are initialized by solving the system of nonlinear equations for a lumped parameter model composed of the Windkessel models used in the multiscale model (figure 3*b*). Thirdly, the solution determined for the lumped parameter model is adapted to compensate for the haemodynamic properties (resistance and compliance) of the multiscale model. As a result, the calibration procedures required a reduced number of calibration iterations. Hence, we conclude that the proposed framework is well suited for a clinical decision setting, where short runtimes are crucial. If required, the computational time can be further reduced by employing modern graphics processing units for massively parallelized computations with a potential speed-up of approximately a factor of 4 [32].

As mentioned above, a major contribution of the framework is that spatially varying personalized regional mechanical wall properties are derived. This is in contrast with previous approaches, where a single PWV value was estimated for the entire domain. As can be observed in figure 5*b*, the spatial variation of the vascular wall properties can be quite pronounced. There are several factors which may contribute to this finding. First of all, all 15 patient datasets of this cohort had aortic valve disease, which leads to a modification of the flow jet through the aortic valve and potentially to an adaptation of the wall properties (especially in the ascending aorta and in the aortic arch). Secondly, previous research has shown that the surrounding tissue of the aorta has a considerable effect on the haemodynamics, leading to higher wave speed, and lower total compliance [33]. Since the various structures surrounding the aorta have different elastic properties, they provide spatially varying external tissue support along the centreline of the aorta.

The main difference from previously introduced methods estimating local vascular wall properties is that herein these wall properties are personalized in the context of a blood flow model, which enforces the physics of fluid flow in elastic domains. Consequently, the results are less affected by measurement noise.

Arterial distensibility is an important factor for the development and assessment of cardiovascular diseases, as elevated systemic vascular stiffness is associated with an increased risk of cardiovascular disease. There are several pathological conditions where the knowledge of regional instead of global mechanical wall properties is of interest. For example, in the case of aortic valve stenosis, a more focused flow jet may affect the aortic wall and cause dilatation of the ascending aorta. Other potential applications are related to the understanding of the development of aortic dissections or aortic aneurysms. Furthermore, previous studies have shown that in the case of aortic coarctation the local compliance plays a crucial role in the estimation of the peak-to-peak pressure drop, a measure that is routinely used for clinical decision making [34]. Furthermore, for stented aortic coarctations, the interplay between the different mechanical properties of the stent and of the surrounding aortic segments is of interest.

The study has a series of limitations. Firstly, it has been evaluated only on 15 patient datasets thus far, and hence the results are preliminary and warrant an evaluation study on a larger number of datasets to be clinically relevant. Secondly, the proposed method has not been tested for anatomical models representative of pathologies which induce severe modifications of the aortic geometry, e.g. aortic coarctation. To account for the effect of the aortic coarctation on the haemodynamics a modified reduced-order model would be required [17]. Thirdly, the reduced-order model introduces an approximation of the geometry, because an axisymmetric, tapering geometry is being considered. However, it has been shown that the one-dimensional model is able to predict time-varying pressure and flow rate waveforms if the tapering is moderate [22]. Fourthly, although we use information about the variation of the cross-sectional areas at many locations of the aorta in our framework, we did not fully take into account global and other movements of the aorta during a heart cycle. Finally, we also note that arterial wall properties are best described using a nonlinear hyperelastic model [35]; however, because the available *in vivo* measurements do not allow for an accurate personalization of such detailed models, we have resorted to a simpler elastic wall model.

## Data accessibility

This article has no additional data.

## Authors' contributions

L.I., D.N., V.M., M.Kr., M.Ke., T.K. and P.S. conceived of and designed the research; L.I., D.N., M.Kr., F.M., M.Ke. and T.K. performed experiments; L.I., D.N., V.M., D.M., M.Kr., M.G. and P.S. analysed data; L.I., D.N., V.M., F.M., M.G. and P.S. interpreted results of experiments; L.I. and F.M. prepared figures; L.I., D.N., V.M. and F.M. drafted the manuscript; L.I., D.N., V.M., F.M. and P.S. edited and revised the manuscript. All authors gave final approval for publication.

## Competing interests

L.I. is an employee of Siemens SRL, Corporate Technology, Brasov, Romania. D.N., F.M. and M.Kr. are employees of Siemens Healthcare GmbH, Erlangen, Germany. V.M., M.G. and P.S. are employees of Siemens Medical Solutions USA, Princeton, NJ, USA.

## Funding

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 611232. This work was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CNCS/CCCDI—UEFISCDI, project number PN-III-P2-2.1-PED-2016-0230, within PNCDI III. Cardioproof is partially funded by the European Commission under the ICT Programme (grant agreement: 611232).

## Disclaimer

This feature is based on research, and is not commercially available. Owing to regulatory reasons its future availability cannot be guaranteed.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3928315.

One contribution of 9 to a theme issue ‘The virtual physiological human: translating the VPH to the clinic’.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.