An optimal control approach to determine resistance-type boundary conditions from in-vivo data for cardiovascular simulations

04/27/2021
by   Elisa Fevola, et al.
0

The choice of appropriate boundary conditions is a fundamental step in computational fluid dynamics (CFD) simulations of the cardiovascular system. Boundary conditions, in fact, highly affect the computed pressure and flow rates, and consequently haemodynamic indicators such as wall shear stress, which are of clinical interest. Devising automated procedures for the selection of boundary conditions is vital to achieve repeatable simulations. However, the most common techniques do not automatically assimilate patient-specific data, relying instead on expensive and time-consuming manual tuning procedures. In this work, we propose a technique for the automated estimation of outlet boundary conditions based on optimal control. The values of resistive boundary conditions are set as control variables and optimized to match available patient-specific data. Experimental results on four aortic arches demonstrate that the proposed framework can assimilate 4D-Flow MRI data more accurately than two other common techniques based on Murray's law and Ohm's law.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 13

page 14

page 15

page 16

10/15/2019

Numerical multiscale methods and effective boundary conditions

We develop numerical multiscale methods for viscous boundary layer flow....
11/04/2019

Reduced order methods for parametric optimal flow control in coronary bypass grafts, towards patient-specific data assimilation

Coronary artery bypass grafts (CABG) surgery is an invasive procedure pe...
08/25/2020

A regularized shallow-water waves system with slip-wall boundary conditions in a basin: Theory and numerical analysis

The simulation of long, nonlinear dispersive waves in bounded domains us...
05/29/2018

Comparison of 1D and 3D Models for the Estimation of Fractional Flow Reserve

In this work we propose to validate the predictive capabilities of one-d...
09/21/2018

Towards a Mini-App for Smoothed Particle Hydrodynamics at Exascale

The smoothed particle hydrodynamics (SPH) technique is a purely Lagrangi...
08/18/2019

A linear system for pipe flow stability analysis allowing for boundary condition modifications

An accurate system to study the stability of pipe flow that ensures regu...
02/24/2021

Learning Off-By-One Mistakes: An Empirical Study

Mistakes in binary conditions are a source of error in many software sys...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction and motivation

In recent years, large improvements have been made in developing patient-specific computational models to predict blood flow in patients affected by various cardiovascular diseases [51, 18]. The availability of data obtained from non-invasive medical imaging techniques, such as computed tomography (CT) and 4D Flow magnetic resonance imaging (MRI), combined with computational fluid dynamics (CFD), has led to more realistic blood flow modeling, with the possibility of investigating the origin and mechanisms behind different cardiovascular diseases.

Mathematical models usually describe blood flow through Navier-Stokes equations, which for complex geometries, such as 3D models of the cardiovascular system, need to be solved numerically. Since it is unfeasible to discretize the entire cardiovascular system, Navier-Stokes equations are solved only on a portion of it, while the rest of the vasculature is accounted for by proper boundary conditions (BCs). Boundary conditions, in fact, model the upstream and downstream vasculatures not included in the 3D model by specifying the physiological conditions at the inlets and outlets of the computational domain of interest. An example of computational model for a human aortic arch can be found in Fig. 1, where the boundary conditions have to be imposed at the inlet and at the four outlets comprising the descending aorta (DAo) and the supra-aortic branches: brachiocephalic artery (BCA), left common carotid artery (LCC), and left subclavian artery (LSUB).

Figure 1. General configuration for the proposed optimal control approach on a patient’s aortic arch. Patient-specific aortic arch haemodynamics is simulated by prescribing the measured inlet flow rate and the optimized resistive-like boundary conditions (BCs) at the outlets, being the outlets the brachiocephalic artery (BCA), the left common carotid artery (LCA), the left subclavian artery (LSUB) and the descending aorta (DAo). These BCs are calculated using pressure and flow measurements at several locations.

The selection of appropriate boundary conditions is always a crucial step when dealing with patient-specific cardiovascular models. Several studies in fact have shown how the choice of boundary conditions has a significant influence on predicted flow patterns as well as on clinically relevant parameters, such as wall shear stress (WSS)[32, 38, 58]. This selection process often benefits from the availability of patient-specific in vivo measurements, so that BCs can be properly tuned to match available clinical data. These data might have been obtained with Phase Contrast Magnetic Resonance Imaging (PC-MRI), a technique for measuring blood velocity in specific sections of the vessel of interest, or with 4D-Flow MRI, a time-resolved PC-MRI with velocity encoding in all three spatial directions and 3D anatomic coverage [49]. For what concerns the inlet boundary condition, usually a velocity profile with a pulsatile waveform is prescribed. If patient-specific data are not available, physiological inlet flow waveforms are taken from the literature [36]. For what concerns outlet BCs, different strategies have been proposed, depending both on the complexity of the treated geometry and on the type of available clinical data. A common choice consists in using constant pressure or traction boundary conditions, even if this greatly affects accuracy [38]. A better solution resorts to the use of zero-dimensional (lumped parameter) models[44, 43, 52, 17, 19], which, by prescribing a specific pressure-flow relationship at each outlet, lead to more physiological pressure and flow waveforms, and to a better approximation of in-vivo data. Lumped models may consist of single resistances or more sophisticated zero-dimensional models, such as the three-element Windkessel model [61]. The use of such outlet boundary conditions, however, leaves the user with the problem of finding appropriate numerical values for the lumped elements. In this sense patient-specific in-vivo measurements, if available, can be used to guide the choice of the parameter values.

The standard approach for cardiovascular CFD simulations consists in a first estimation of outlet resistive BCs through Murray’s law [33], which gives an indication of how blood flow splits in vessel bifurcations. These initial estimates are then refined to match available patient-specific data through a manual tuning process, until the user-desired accuracy is obtained [5, 31]. Despite being straightforward to implement, this approach is sub-optimal, operator-dependent, time-consuming, and affected by non-repeatability. To achieve robust and reliable patient-specific simulations, it is vital to rely on automated procedures for the selection of parameter values [6, 48]

. A first class of solutions is based on the use of a Kalman filter and its non-linear extensions, which performs the assimilation of measurements

on the fly and updates the unknown parameters accordingly [27, 7]. Arthurs et al. [4] proposed a reduced-order unscented Kalman filter for the sequential estimation of simple lumped parameter network parameters. Since the uncertainty estimation was limited to the unknown parameters, and not to the entire state space, the filtering operation was computationally tractable directly on 3D models. The data assimilation process, however, is still computationally expensive, as the estimation of n parameters requires running n+1 forward simulations.

A different approach for parameter estimation relies on Bayesian inference, where the unknown parameters are described as random variables, with the goal of yielding confidence regions instead of point estimates and quantifying the uncertainty affecting simulation results 

[56, 45, 13, 37]. Schiavazzi et al. [45], for example, used Bayesian estimation for automated tuning of 0D lumped model parameters to match clinical data, estimating also the uncertainty induced by errors in the measurements. Tran et al. [56]

proposed an iterative approach based on adaptive Markov chain Monte Carlo sampling to quantify confidence in local hemodynamic indicators, e.g., wall shear stress and oscillatory shear index (OSI).

Alternatively, variational approaches, such as optimal control, perform the parameter estimation by treating the governing equations as state system and the mismatch between the state solution and the patient-specific measurements as cost functional to be minimized [25, 62, 16, 12]. In practice, the unknown boundary conditions are chosen as control variables to minimize the cost functional. As other advanced data assimilation methods, also optimal control suffers from a high computational cost, since it requires solving a constrained minimization problem. Nevertheless, recent studies have documented its successful application to several haemodynamics problems [40, 10, 34]. For what concerns the calibration of hard-to-quantify boundary conditions, Tiago et al.[54], for example, proposed a velocity control approach for the assimilation of velocity measurements in blood flow simulations. The framework was validated on an idealized 2D geometry and tested on a 3D model of a brain aneurysm, but with the assimilated velocity data still synthetically generated. In the work by Koltukluouglu et al. [25], a data assimilation method based on optimal boundary control for 3D steady-state blood flow simulations was presented. The methodology was validated against real 4D-flow MRI data measured on a glass replica of a human aorta. To reduce the computational cost of optimal control, Romarowski et al. [42] identified a surrogate optimization problem for prescribing PC-MRI data as outlet boundary conditions, tuning the parameters of a three-element Windkessel model via a least-squares approach. Zainib et al. [62] proposed a reduced order framework for the application of optimal control to coronary artery bypass grafts, where synthetically-generated measurements were used to tune Neumann-type outlet boundary conditions, thus bringing optimal control closer to a real clinical setting. Tiago et al. [53], moreover, used an optimal control approach to address the uncertainty coming from the segmentation process, which affects the definition of the geometry and, in turn, wall shear stress and its derived measures.

In this work, we propose an automated method for the selection of resistive-type outlet boundary conditions based on the solution of an optimal control problem, which assimilates a set of patient-specific measurements, e.g., through 4D-Flow MRI imaging. The assimilation is performed by constraining the minimization of the objective functional to the solution of a steady Stokes problem. Since Stokes equations correspond to the linearized steady-state Navier-Stokes equations, they are less computationally expensive, while still providing an adequate estimation of boundary conditions. The control variables are the unknown resistances imposed at the outlet through the coupled multidomain method introduced by Vignon-Clementel et al.[59]. To ensure a better match between simulation results and in-vivo data, the inlet flow waveform measured with 4D-Flow MRI is imposed as a patient-specific inlet boundary condition, as reported in Fig. 1. One of the novel aspects of this work is the use of more realistic boundary conditions with respect to Neumann and Dirichlet BCs, up to now the standard in an optimal control setting. Its suitability for real clinical scenarios is demonstrated in this work by validating the framework on four aortic arches, reconstructed from medical images of real clinical cases. The presented method is used in these cases to set the outlet BCs on the descending aorta and the supra-aortic branches, assimilating real 4D-Flow MRI data. The boundary conditions obtained with the proposed method are compared to those obtained with two alternative calibration techniques, namely, Murray’s law and Ohm’s law, demonstrating the ability of optimal control to assimilate known physiological data consistently better. Moreover, an analysis of time-averaged wall shear stress and oscillatory shear index values obtained using the three different calibration methods is used to assess their effect on clinical haemodynamic indicators.

This work is organized as follows: Section 2 presents the methodological details of the proposed optimal control approach. Results on four patient-specific anatomies are reported in Section 3, while Section 4 analyses the advantages and limitations of the proposed approach and Section 5 provides our conclusions and future perspectives.

2. Methodology

This section presents the methodological aspects of the proposed framework based on optimal control. Before solving the actual optimal control problem, some preliminary steps are required, namely, the extraction of patient-specific cardiovascular configurations from clinically acquired images, and the acquisition of patient-specific 4D-Flow MRI data. The proposed framework is represented in Fig. 2.

Figure 2. Scheme of the proposed framework for boundary condition estimation through optimal control.

2.1. Data acquisition and anatomical reconstruction

The proposed methodology has been tested in four clinical cases from a single-center prospective study conducted at the Sunnybrook Health Sciences Centre in Toronto, Canada. The study was approved by the local ethics board and informed consent was obtained. Moreover, all the measurements were acquired non-invasively. Patients presented at the hospital for coronary bypass graft surgery. Between three and six weeks after surgery, a cardiac computer tomography (CT) was performed and anatomical information about their aorta and its supra-aortic branches was acquired using 320-detector row CT scanner (Aquilion One, Canon Medical Systems). From CT images, the vessels surface was reconstructed using the open-source package SimVascular 

[57]. The reconstructed volume was discretised into tetrahedral elements using TetGen [46]. After the CT scan, the blood velocity in the aorta and its branches was acquired in-vivo using a 4D-flow magnetic resonance imaging (MRI) sequence, using a 3T MRI scanner (MAGNETOM Prisma, Siemens Healthineers). The acquisition was performed using a 4D flow imaging sequence with retro-gating and adaptive navigator respiratory gating. Imaging parameters were as follows: velocity encoding=150 cm/s, field of view=200-420 mm x 248-368 mm, spatial resolution=1.9-3.5 x 2.0-3.2 x 1.8-3.5 mm, temporal resolution=39.9-47.2 ms, flip angle=8. After 4D-Flow MRI, diastolic blood pressure and systolic blood pressure were measured with brachial cuff-based method. The mean arterial pressure was computed as

(1)

2.2. Determination of boundary conditions through optimal control

In general, an optimal control problem consists in determining a number of control variables, which are unknown quantities, so that a cost functional is minimized, under the constraint of a system of equations describing the physical state of blood flow. The solution of this problem provides both an estimate for the unknown control variables, and the solution of the state variables (pressure and velocity). In the proposed framework, blood flow is modeled with the Stokes equations. This choice ensures a simpler optimal control problem with respect to its nonlinear version based on Navier-Stokes equations. Given the simplification that the Stokes equations introduce on the model, we will assess the effect on the quality of the determined output boundary conditions, by validating the obtained results on a Navier-Stokes model. This analysis is reported in Section 3.

Considering the large diameter of the aorta and supra-aortic branches, we can safely assume that the blood behaves as a Newtonian fluid and that the viscosity is constant, since the dimension of blood particles is smaller compared to vessel diameters. Therefore, we adopt the incompressible Stokes equations as state equations modeling blood flow in the computational domain , which is the volume of the 3D aortic arch model reconstructed from medical images, as shown in Fig. 1. We refer to the boundaries of as , where , , and denote the inlet of the aorta, the vessel walls, and the outlets of the aortic arch, with ( in the case of Fig. 1), respectively. The control variables are denoted as . Thus, state equations can be written in strong form as

(2)

where is blood velocity, the pressure, dynes/cms the dynamic viscosity, and the outward normal to the outlets. A plug profile is imposed at the inlet whose average value is extracted from the 4D-Flow MRI data. A no-slip condition is imposed at the vessel walls , assumed to be rigid and non-permeable, while at the outlets a resistive type boundary condition is imposed, following the coupled multidomain method proposed by Vignon-Clementel et al [59]

. Since the optimal control framework for partial differential equations 

[21] is typically derived starting from the weak formulation of the state equations, we introduce Hilbert spaces and for the velocity and the pressure , respectively. In particular, we choose and . Furthermore, we denote by the space associated to the controls . We note that and are function spaces (i.e., and are functions depending on the spatial coordinate ), while is simply an Euclidean space (i.e., are scalar numbers, and not functions).

Starting from the strong form in (2), the weak formulation is derived as: given , find

such that

(3)

for every

where and are the test functions associated to velocity and pressure, respectively.

We now reformulate the term

in to obtain an equivalent weak formulation which is more suitable for the forthcoming finite element discretization. We introduce a set of Lagrange multipliers , , defined as

We further denote by

the vector collecting the Lagrange multipliers, and

its associated space. Therefore, the equivalent weak formulation is: given , find such that

(4)

for every , where collects the test “functions” associated to the Lagrange multipliers . This is the final system representing how the control affects the underlying physics of the model. In order to set up the optimal control problem, a proper cost functional must also be defined.

We define the cost functional with the following form

(5)

The first term in (5) represents the normalized difference between the state pressure and the patient’s average pressure measured at a specific cross section ; as explained in 2.1, in this work was assumed equal to the mean arterial pressure, computed from the measured systolic and diastolic pressure. The second term, instead, represents the normalized difference between the calculated flow rate (obtained integrating velocity on the outlet section) and the flow rate extracted from 4D-Flow MRI data at each outlet . The choice of minimizing the normalized difference between simulated and measured quantities ensures that each term gives the same contribution to the optimization process, even when assimilating measurements with different orders of magnitude. The weights and , however, can be used to change the contribution of each term individually. For the aortic arches under analysis, and all were set to 1, meaning that all measurements contribute equally to the minimization process. The optimal control problem then reads:

Problem 1.

Find such that functional (5) is minimized, under the constraint that satisfy Equation (4).

To solve the optimal control problem, we adopt the adjoint-based Lagrangian approach [39, 22, 21, 23]. This approach models the optimal control problem as an unconstrained minimization problem, whose solution corresponds to the minimum of a properly defined Lagrangian functional. In practice, the optimal solution is the one where all the derivatives of the Lagrangian functional vanish. The Lagrangian formulation requires the introduction of three new unknowns, the so-called adjoint variables, , and . The Lagrangian functional for this problem takes the form

(6)

Given (6), the optimality system can be obtained by imposing that the derivatives of with respect to must vanish, in short requiring . Taking the derivative of (6) with respect to (denoted by ) in the direction we obtain

(7a)
while taking the derivative of (6) with respect to in the direction we get
(7b)
Similarly, the derivatives of (6) with respect to are, respectively,
(7c)
(7d)
(7e)
(7f)
(7g)

The presence in Eq. (7a) of a term containing the product of two integrals requires the use of an additional Lagrange multiplier. We thus introduce the new variables , which we substitute in Eq. (7a), and we add the following equations to the system

(8)

Equations (7a) through (8) form the so-called coupled optimality system, which we solve through a one-shot approach [21, 50], where the system is solved directly for all the unknown variables.

We adopt an optimise-then-discretise approach, meaning that we introduce the numerical discretisation of the problem after the derivation of the optimality system. In particular, we rely on Galerkin finite element methods for the solution of the discrete version of the system. The domain was discretised into a finite mesh of size , and so we introduce finite-dimensional solution spaces . In particular, we use Taylor-Hood elements for the velocity-pressure pair, i.e. finite elements to define and finite elements for . Since the spaces associated to control and Lagrange multipliers are already finite-dimensional, we set and . The discretised system is solved using the open-source libraries FEniCS [2, 29] and multiphenics [1], the latter being an open-source library developed at SISSA mathLab for easy prototyping of problems characterized by a block structure and boundary restricted variables. The numerical solution of the problem is obtained by means of MUMPS [3], a parallel sparse direct solver.

3. Numerical Results

3.1. Synthetic data

A mesh convergence analysis was first conducted on the patient-specific meshes used for the experiments. In Figure 5 the average TAWSS and OSI are plotted with respect to the number of elements in the mesh. The mesh used in the following experiments has about elements, which has reached convergence for OSI, but not yet for WSS. This choice, however, is the best compromise in terms of accuracy and computational cost of the simulations.

The validity of the proposed approach was first tested on a test case (corresponding to the patient anatomy that will be labeled as “case 1” in the following) with synthetically generated data. To obtain the ground truth data, case 1 was simulated in SimVascular, imposing a velocity waveform at the inlet with average flow rate , using a blunt profile. At the outlets, three-elements Windkessel models were imposed. A physiological value of the total resistance at each outlet was arbitrarily chosen, and each resistance was then split into a proximal one, , and a distal one, , as suggested by Kim et al [24]. For the capacitance, a total value of cm/dyn was assumed [28], which was then split among the four outlets proportionally to their area. The distal pressure was set to 0 at all outlets. To reach periodic convergence, the simulation was run for five cardiac cycles, and the average flow rate at the four outlets and the average pressure were extracted from the last cardiac cycle. These data were fed into the optimal control tool presented in Section 2 to estimate the total outlet resistances. The estimated parameters were then used to set the outlet boundary conditions of a second simulation, again in Simvascular. Both resistances and capacitances were split adopting the same rules of the forward simulation. Table 1 reports a comparison between original and estimated resistance values, together with original and estimated average flow rates at the outlets. Flow and pressure waveforms at the outlets are compared in Fig. 10, which shows how the resistances chosen with optimal control allow to reconstruct the original flow waveform with a relative error of 0.09% for BCA, 0.09% for LCC, 0.1% for LSUB and 0.02% for DAo. The pressure waveform is recovered with a relative error of 0.005%.

Resistance (dyns/cm) Flow rate (cm/s)
BCA LCC LSUB DAo BCA LCC LSUB DAo
Ground truth 7,000 21,000 16,000 1,700 19.02 6.30 8.13 80.02
Estimated 6,937 20,846 15,991 1,685 19.05 6.30 8.09 80.03
Table 1. Results of experiment with synthetic data.
(a) Convergence analysis of mean TAWSS.
(b) Convergence analysis of mean OSI.
Figure 5. Mesh convergence analysis.
(a) Flow rate waveform at Brachiocephalic artery (BCA)
(b) Flow rate waveform at Left Common Carotid artery (LCC)
(c) Flow rate waveform Left Subclavian artery (LSUB)
(d) Pressure waveform at Descending Aorta (DAo)
Figure 10. Comparison of flow rate waveforms at the three supra-aortic branches of case 1, and of pressure waveform at the Descending Aorta, for the test case with synthetically generated data.

3.2. Patient-specific data

Moving to the real cases with patient-specific data, the outlet boundary conditions estimated with the proposed optimal control framework were compared to those obtained with two alternative common techniques based on Murray’s law and Ohm’s law.

Murray’s law

Murray’s law [33], formulated by Cecil D. Murray in 1926, governs the branching pattern of vessels, such that the flow in each outlet is proportional to its cross-sectional area. In particular, the general form of Murray’s law reads

(9)

where is the flow rate at the outlet, is the radius at the outlet, and for the aortic arch is conventionally set to 2. Multiplying by to express the relationship in terms of outlet areas, and using , one can estimate the outlet resistances for the aorta and its main branches as

(10)

where is the sum of the area of all the aortic outlets, while is the area of the outlet to which resistance is associated. The total resistance was computed as the mean pressure measured non-invasively on the patient, as reported in Section 2.1, divided by the mean aortic flow rate measured with 4D-Flow MRI, and then split among the outlet branches according to Equation (10). The application of Murray’s law for predicting flow splitting has been largely investigated both on human and animal subjects by a number of studies [20, 55, 8, 41], which evidenced its validity on a large portion of the cardiovascular system, even if with some limitations on the first branches of the aortic arch [63]. Note that, except for the average pressure and inlet flow rate, Murray’s law does not take into consideration patient-specific measurements to compute flow splitting: instead, assuming that branches with larger cross-section have higher flow rates, it is solely based on the patient’s anatomy. This feature justifies its use in those studies where in-vivo measurements of flow rates are not available [47, 60, 15].

Ohm’s law

Based on the Ohm’s law, it is possible to set outlet resistances by taking advantage of the analogy between the cardiovascular system and electrical circuits. In particular, knowing the mean pressure computed from diastolic and systolic pressure measured non-invasively and the outlet flow rates measured with 4D-Flow MRI, outlet resistances can be computed as

(11)

This method is also common [38] and it is based on the idea of performing the parameter estimation on 0D models [35] but, differently from Murray’s law, it requires the availability of outlet flow rates measured in-vivo. The measured flow rates usually respect the mass conservation principle with approximately a deviation [11], due to measurement uncertainty. For this reason, we developed a minimization problem which, using Ohm’s law, estimates outlet resistances while trying to impose the mass conservation principle. In this way, we try to compensate for the intrinsic inconsistency in the data, which is not accounted for in Ohm’s law of Eq. (11). This is achieved by approximating the aorta with its equivalent 0D model, which is represented in Fig. (a)a, and minimizing the cost function

(12)

The cost function was minimized using the Matlab function fminsearch, which employs a derivative-free simplex search method [26]. The cost functional reported in Equation (12) deliberately replicates the one used inside the optimal control framework, reported in Equation (5), with the first term representing the normalized difference between the computed pressure and the measured one (), and the second the outlet flow rates . However, optimal control relies on Stokes equations as a 3D-model of the underlying system, whereas here a 0D approximation is used.

(a)

(b)
Figure 13. Adopted equivalent circuits. On the left, equivalent circuit used for Ohm’s law method. On the right, 3-element Windkessel model used as outlet boundary condition for unsteady Navier-Stokes simulations.
Figure 14. Velocity streamlines (top row) and pressure distributions (bottom row) for case 3 at three different time instants. The time points T1 (mid-systolic acceleration, left column), T2 (peak systole, middle column) and T3 (diastole, right column) are defined along a flow waveform shown on the left.

The experiments described in this section were conducted on four patient anatomies, obtained as described in Section 2.1. Table 2 reports for each case the measured flow rates. Specifically, the column labeled Net flow indicates the difference between measured inlet flow and outlet flows, thus quantifies the violation of mass conservation on the measurements. The inlet flow reported in Table 2 was obtained by subtracting 4% to the flow measured in the ascending aorta, which corresponds to the total coronary circulation [14], not considered in the models. The presence of flow rate inconsistencies in 4D-Flow MRI data could be due to the finite resolution of 4D-flow MRI, motion artifacts, and the presence of noise, especially in presence of complex helical and vortical flows [9]. The amount of net flow is below for all the cases under analysis, which is considered acceptable for 4D-Flow MRI measurements according to Dyverfeldt et al [11]. While in case 1 and case 4 the net flow is significant, for cases 2 and 3 it is practically negligible.

3.3. Boundary conditions estimation: Stokes model

The main results are reported in Table 3 where, for each case, the first row reports the patient’s average pressure, measured non-invasively after MRI, and the flow rates measured in-vivo with 4D-Flow MRI at the four outlets: BCA, LCC, LSUB and DAo. The values obtained solving Stokes equations on the 3D geometry, with resistance values calculated by means of Ohm’s law (Eq. (11)) are reported in the second row, while the third row contains the results obtained with the optimization based on Ohm’s law as in Eq. (12). The fourth contains the results for Murray’s Law. Finally, the fifth row shows the results of the proposed methodology based on optimal control. In all simulations, the inlet flow was imposed using a plug profile. For case 2 and case 3, where the net flow is negligible, both Ohm’s law and optimal control properly assimilate the available data, with a relative error of less than on all the outlet flow rates. For case 1 and case 4, where the net flow is high, optimal control outperforms the other techniques, achieving the smallest relative errors on pressure and BCA, LCC, LSUB flow rates. In presence of a large net flow, we cannot expect a perfect assimilation of measurements, as they are intrinsically non-physical. In that case, optimal control shows a smaller sensitivity to inconsistencies in the data, leading to a physical solution which is closer to measurements with respect to the other techniques. As expected, the optimization based on Ohm’s law is more accurate than Ohm’s law for case 1 and 4, while the two methods are basically equivalent for case 2 and 3. As already pointed out earlier, as Murray’s law is the only technique which does not take into account the measured flow rates to set outlet boundary conditions, the corresponding solution is the one that most deviates from patient measurements. The solution of the optimal control problem required an average of 6.75 minutes (wall clock CPU time), running on 18 Lenovo SD530 nodes, each with 40 Intel ”Skylake” cores and 202 GB RAM.

Case
number
Inlet flow
(cm/s)
Outlet flow
(cm/s)
Net flow
(cm/s)
Mass conservation
violation (%)
1 119.10 103.46 15.64 13%
2 107.00 107.18 -0.18 0.17%
3 125.63 125.60 0.03 0.02%
4 103.00 90.21 12.79 12.4 %
Table 2. Table summarizing, for each case, the measured inlet flow rate, the sum of the measured outlet flow rates, the difference between the two (Net flow), and the percentage of mass conservation violation.

max width = Case number Method Pressure (mmHg) Flow rate (cm/s) (% error w.r.t. measurements) BCA LCC LSUB DAo 5 Measurements 98.7 15.9 5.98 8.48 73.1 Ohm’s law 113 (15%) 18.3 (15.3%) 6.89 (15%) 9.77 (15%) 84.2 (15%) Opt. based on Ohm’s law 108 (9.5%) 17.6 (11%) 6.55 (9.5%) 9.28 (9.4%) 85.8 (17%) Murray’s law 98.8 (0.1%) 19.2 (20%) 6.19 (3.5%) 7.48 (-11%) 86.2 (18%) Proposed 98.7 (0%) 16.6 (4.4%) 6.08 (1.7%) 8.67 (2.2%) 87.8 (20%) 9 Measurements 105 13.2 6.71 7.47 79.8 Ohm’s law 105 (0%) 13.1 (-0.7%) 6.68 (-0.4%) 7.44 (-0.4%) 79.5 (-0.4%) Opt. based on Ohm’s law 105 (0%) 13.1 (-0.7%) 6.70 (-0.14%) 7.45 (-0.26%) 79.6 (-0.3%) Murray’s law 106 (0.9%) 17.0 (29%) 6.28 (6%) 9.30 (24%) 74.2 (-7%) Proposed 106 (0.9%) 13.2 (0%) 6.71 (0%) 7.47 (0%) 79.5 (-0.4%) 10 Measurements 103 19.0 11.3 10.5 84.8 Ohm’s Law 103 (0%) 19.0 (0%) 11.3 (0%) 10.5 (0%) 84.9 (0.1%) Opt. based on Ohm’s law 104 (0.65%) 19.1 (0.5%) 11.4 (0.62%) 10.5 (0%) 84.9 (0.1%) Murray’s law 104 (0.65%) 18.1 (-5%) 9.37 (-17%) 7.22 (-31%) 91.2 (7%) Proposed 103 (0%) 19.0 (0%) 11.3 (0%) 10.5 (0%) 85.0 (0.2%) 11 Measurements 100 9.87 4.32 6.92 69.1 Ohm’s law 115 (15%) 11.3 (14%) 4.94 (14%) 7.91 (14%) 78.9 (14%) Opt. based on Ohm’s law 109 (9%) 10.8 (9%) 4.68 (8.3%) 7.49 (8.2%) 80.1 (16%) Murray’s law 101 (1%) 15.3 (55%) 4.09 (-5%) 5.09 (-26%) 78.6 (13.7%) Proposed 100 (0%) 10.1 (2.6%) 4.37 (1.1%) 7.04 (1.7%) 81.5 (18%)

Table 3. Comparison of pressure and flow rates for Stokes simulations with boundary conditions obtained with the different methods.
Case
number
Method Resistance (dyns/cm)
BCA LCC LSUB DAo
1 Ohm’s law 8,288 21,979 15,518 1,800
Opt. based on Ohm’s law 8,190 21,981 15,511 1,679
Murray’s law 6,837 21,242 17,591 1,527
Proposed 7,941 21,609 15,153 1,497
2 Ohm’s law 10,690 21,003 18,847 1,764
Opt. based on Ohm’s law 10,737 20,987 18,867 1,767
Murray’s law 8,306 22,451 15,161 1,900
Proposed 10,699 21,006 18,864 1,773
3 Ohm’s law 7,248 12,142 13,094 1,624
Opt. based on Ohm’s law 7,255 12,145 13,090 1,629
Murray’s law 7,646 14,723 19,102 1,513
Proposed 7,249 12,131 13,078 1,621
4 Ohm’s law 13,600 31,060 19,391 1,943
Opt. based on Ohm’s law 13,500 31,069 19,399 1,815
Murray’s law 8,751 32,800 26,358 1,708
Proposed 13,166 30,500 18,903 1,634
Table 4. Resistance values chosen by the different methods.

3.4. Boundary conditions estimation for inaccurate measurements

In presence of a large net flow, the presented approach compensates for the inaccuracy in the data by adjusting the outlet boundary conditions. This means that the inconsistencies in the measurements are resolved entirely at the outlets, while the measurement imposed at the inlet is considered deterministic. It would be desirable, instead, to gain some flexibility in the assimilation of the inlet flow, which is equally affected by uncertainty. In this case, a modification of the approach presented in Section 2 is possible, which estimates both the inlet and outlet BCs by means of an additional control at the inlet Dirichlet boundary condition. In particular, the velocity at the inlet in 2 is expressed as

(13)

where is a scalar control variable and is the inlet velocity profile, with average value equal to the one measured with 4D-Flow MRI. By choosing the best value for , the optimal control problem will be able to change the inlet flow rate to better assimilate the available data. This obviously requires a slight modification of the cost functional, with an additional term for assimilating the flow rate at the inlet :

(14)

Moreover, the Dirichlet control requires weakly imposing the inlet condition by means of Lagrange multiplier, thus increasing the final size of the system of equations (7a) - (8). In table 5 we report the results with this alternative formulation for cases 1 and 4, which had the largest Net flow. As expected, the net flow is now resolved acting on all boundary conditions, leading to lower differences between measured and simulated flow rates. It is worth mentioning that the additional control leads to an increase in the dimensions of the problem, and consequently to larger computational costs (an average of 10 minutes of wall clock CPU time, running on 24 Lenovo SD530 nodes, each with 40 Intel ”Skylake” cores and 202 GB RAM).

(a) case 1
(b) case 2
(c) case 3
(d) case 4
Figure 19. Comparison of pressure and outlet flow rates for time-dependent Navier-Stokes simulations. The histograms report the relative difference with respect to the corresponding patient measurements.

max width = Case number Method Pressure (mmHg) Flow rate (cm/s) (% error w.r.t. measurements) Inlet BCA LCC LSUB DAo 1 Measurements 98.7 119.1 15.9 5.98 8.48 73.1 Proposed 98.7 (0%) 107.9 (-9%) 16.1 (1.2%) 6.01 (0.5%) 8.53 (0.6%) 77.3 (5.7%) 4 Measurements 100 103.0 9.87 4.32 6.92 69.1 Proposed 100 (0%) 94.3 (-8%) 9.95 (0.8%) 4.33 (0.23%) 6.96 (0.6%) 73.0 (6%)

Table 5. Comparison of measurements and pressure and flow rates for Stokes simulations obtained with the extended optimal control formulation, which controls both the inlet and the outlets boundary conditions.

3.5. Unsteady Navier-Stokes simulations

An accurate representation of blood flow in the aorta is generally obtained through unsteady Navier-Stokes simulations, which provide a more realistic and accurate time evolution of blood flow. Assuming laminar regime, the incompressible Navier-Stokes equations take the form

(15)

We used the resistance values reported in Table 4 as outlet boundary conditions of high-fidelity unsteady Navier-Stokes simulations performed using SimVascular [57]. The purpose of this analysis is twofold. First, it verifies that the flow splitting obtained with a steady linear Stokes model is still valid in a non-linear, unsteady scenario. Second, it allows to analyse the impact that boundary conditions obtained with different estimation techniques have on wall shear stress-related indicators, which are clinically relevant. For unsteady Navier-Stokes simulation, the type of outlet boundary condition which best replicates realistic flow conditions is the three-element Windkessel model [38] represented in Fig. (b)b. Referring to Fig. (b)b, each resistance was split into a proximal one, , and a distal one, , as suggested by Kim et al [24]. For the capacitance, a total value of cm/dyn was assumed [28], which was then split among the four outlets proportionally to their area. The distal pressure was set to 0 at all outlets. The dynamic viscosity was set to dynes/cms, a rigid wall model was assumed and a plug profile was imposed at the inlet. The inlet waveform was the one extracted from 4D flow MRI. The time-step value for the transient simulations was set to 0.5 ms, and 5 cardiac cycles were simulated. The results reported here refer to the last cardiac cycle. Each simulation required an average of 8 hours (clock wall CPU time), running on 4 Lenovo SD530 nodes, each with 40 Intel ”Skylake” cores and 202 GB RAM.

Figure 14 reports the velocity streamlines and pressure distribution for case 3 at three different time instants along the cardiac cycle. A comparison of pressure and outlet flow rates obtained with Navier-Stokes simulations using outlet boundary conditions estimated with the three different techniques introduced previously is reported in Fig. 19. The histograms represent, for each outlet, the relative difference of the average flow rate with respect to the measured one using Murray’s law, Ohm’s law, or the proposed method. Also for this simulations, the inlet flow was imposed using a plug profile. As expected, the errors of the obtained flows with respect to the measured ones increased when moving from a Stokes model to a Navier-Stokes one, mostly due to the non-linearity and time-dependency introduced by the latter. Nevertheless, similar trends to those of the Stokes experiments reported in Table 3 can be observed when moving to Navier-Stokes simulations. In particular, Murray’s law remains the method providing the largest deviations from measured flow rates, while both Ohm’s law and optimal control are closer to assimilated data. With the exception of case 3, optimal control is still the method which best replicates measured flow rates. These results show that the BCs estimated with a linear, steady Stokes model proved to be a good choice when moving to high-fidelity Navier-Stokes simulations, thus supporting the approach of estimating BCs on a linearised Stokes model.

Additionally, Figure 24 reports a comparison of time-dependent flow rates extracted from 4D-Flow MRI and those obtained from time-dependent simulations with boundary conditions estimated with optimal control. Even if the proposed estimation method assimilates only the average flow rates, the time-dependent waveforms are still recovered with a good degree of accuracy. For the sake of space, only results for case 4 are reported, with comparable results for the other cases.

(a) Flow rate waveform at Brachiocephalic artery (BCA)
(b) Flow rate waveform at Left Common Carotid artery (LCC)
(c) Flow rate waveform Left Subclavian artery (LSUB)
(d) Flow rate waveform at Descending Aorta (DAo)
Figure 24. Comparison of 4D-Flow MRI flow waveforms and simulated flow waveforms with boundary conditions estimated with the proposed optimal control approach for case 4.

To assess the influence that the adopted BCs estimation technique may have on clinically relevant parameters, we carried out an additional analysis on two relevant haemodynamic indicators, namely, the time-averaged wall shear stress (TAWSS) and the oscillatory shear index (OSI), calculated from Navier-Stokes simulation results using the equations reported by Martin et al. [30]. For the sake of space, we report the analysis for case 3, but similar results were obtained for the other cases. Fig. 25 shows, in the left column, the TAWSS obtained for the two cases with three different techniques (optimal control, Murray’s law, and optimization based on Ohm’s law), and in the right column the local relative difference with respect to optimal control results. The same analysis was repeated for the OSI in Fig. 26. For each point of the surface anatomy, the local relative difference was computed as:

(16)
(17)

For case 3, the difference in the TAWSS reaches a maximum relative difference of , while the discrepancy in the OSI value reaches . The largest differences occur in the Murray’s law case, in correspondence of LSUB, which is also where the estimated resistance values differ the most (). It is worth noticing that, outside some specific ‘hot-spots’, the relative error is generally lower, around 10%. This analysis reveals the impact that the values of resistive boundary conditions have on haemodynamic indicators. In particular, given the relevance of TAWSS and OSI in a clinical context, the adoption of different techniques for BCs estimation could possibly affect the observations done by medical doctors, reaffirming the importance of adopting an automated, reliable, and operator-independent technique for boundary condition estimation.

Figure 25. Left: Time-averaged wall shear stress (TAWSS) results obtained with the three different BC estimation techniques for case 3. Right: local absolute differences in TAWSS as compared to results obtained with optimal control. Anterior and posterior views of the anatomy are provided.
Figure 26. Left: Oscillatory Shear Index (OSI) results obtained with the three different BC estimation techniques for case 3. Right: local absolute differences in OSI as compared to results obtained with optimal control. Anterior and posterior views of the anatomy are provided.

4. Discussion and Limitations

The proposed approach for boundary condition estimation presents some considerable advantages. First, optimal control relies on a rigorous and general mathematical framework, which ensures the possibility to apply it to different vessels in the cardiovascular system, provided that in-vivo measurements are available. Second, the mathematical structure of optimal control leaves great flexibility in terms of the nature and quantity of measurements to assimilate, an appealing feature for cardiovascular applications, where in-vivo measurements can be either pressure or flow waveforms. In particular, the proposed method does not specifically require data from 4D-Flow MRI, but it is compatible with any method providing flow rate information at the inlets and outlets of the region of interest, including the more common PC-MRI modality. Moreover, the proposed formulation assimilates data in a least-squares sense [21], reducing the influence of stochastic measurement uncertainty on the solution.

On the other hand, the proposed technique presents some limitations. The main limitation comes from the use of steady-state, linearized Navier-Stokes equations inside the optimal control problem. As explained in the Introduction, optimal control problems are characterized by a high computational cost. For this reason, the solution of a time-dependent optimal control problem is at the moment computationally intractable. The necessity to use a steady-state formulation of the problem requires resorting to the Stokes model, as Navier-Stokes equations may not converge to a steady-state solution for a complex flow like the one in the aortic arch. On one side, using a simplified model like the Stokes one significantly reduces the computational cost for the data assimilation process. On the other side, steady Stokes equations are clearly inadequate to get a realistic representation of the flow in the aorta. For this reason, the solution of the optimal control problem is used solely to estimate the outlet resistances, while for realistic pressure and velocity distributions a subsequent time-dependent simulation is still necessary. One main consequence of using a steady-state formulation is the possibility to estimate only the total value of the resistance. Estimating the outlet capacitance, in fact, would require solving a time-dependent Navier-Stokes problem. For the same reason, it is not possible to determine multiple resistance values at each outlet when using more advanced outlet boundary conditions, such as the proximal and distal resistances in a 3-element Windkessel model, being then forced to adopt capacitance values and rules for resistance splitting taken from literature. Lastly, the proposed method can only assimilate average flow rates, and not time-dependent flow waveforms. The results presented in Section 

3, in particular the simulated waveforms reported in Figure 24, prove that, when the resistance and capacitance values determined with the proposed method are used in an unsteady simulation, the obtained pressure and flow rates values are in sufficient agreement with the in-vivo measurements used for BC estimation.

5. Summary and Conclusion

In this work, we proposed a framework based on optimal control for the automated estimation of unknown resistance-type boundary conditions, while assimilating in-vivo pressure and flow rate measurements. The experiments conducted on four patient anatomies revealed the validity of the proposed optimal control-based technique for assimilating 4D-MRI data. Specifically, when compared to two other common techniques, namely, Murray’s law and Ohm’s law, the proposed framework performed consistently better. An additional investigation of the effects of the different estimation methods on wall shear stress-related parameters exposed the influence that boundary conditions play on clinically relevant quantities. This further confirms the need for an automated method, such as the one proposed, which eliminates the expensive manual tuning phase, together with intra- and inter-operator variability. The proposed method, therefore, represents a first step in incorporating optimal control, thus a reliable, automated, and robust optimization technique into a framework which can be exploited by the medical community in a clinical setting. The field is promising and opens important perspectives for mathematical modelling and numerical simulation in cardiovascular flows.

A first extension of the present work is the application to other clinically relevant scenarios, such as coronary artery bypass grafts, where the scarcity and noise of available data limit the applicability of other estimation techniques. Furthermore, the current framework estimates only resistive-type quantities, while the capacitance values are still chosen based on generic information available in the literature, usually combined with a manual tuning process. Future works will be directed toward the automated tuning of capacitance values as well.

Acknowledgements

We acknowledge the support provided by the Canada Research Chairs program, the Radiological Society of North America, Compute Canada, and the H2020 ERC CoG AROMA-CFD project (G.A. 681447).

References

  • [1] multiphenics – easy prototyping of multiphysics problems in fenics. https://mathlab.sissa.it/multiphenics, 2021. Accessed: 21.03.2021.
  • [2] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E Rognes, and Garth N Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015.
  • [3] Patrick R Amestoy, Iain S Duff, Jean-Yves L’Excellent, and Jacko Koster. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications, 23(1):15–41, 2001.
  • [4] Christopher J Arthurs, Nan Xiao, Philippe Moireau, Tobias Schaeffter, and C Alberto Figueroa. A flexible framework for sequential estimation of model parameters in computational hemodynamics. Advanced Modeling and Simulation in Engineering Sciences, 7(1):1–37, 2020.
  • [5] Alessandro Boccadifuoco, Alessandro Mariotti, Katia Capellini, Simona Celi, and Maria Vittoria Salvetti. Validation of numerical simulations of thoracic aorta hemodynamics: comparison with in vivo measurements and stochastic sensitivity analysis. Cardiovascular Engineering and Technology, 9(4):688–706, 2018.
  • [6] Mirko Bonfanti, Gaia Franzetti, Gabriele Maritati, Shervanthi Homer-Vanniasinkam, Stavroula Balabani, and Vanessa Díaz-Zuccarini. Patient-specific haemodynamic simulations of complex aortic dissections informed by commonly available clinical datasets. Medical Engineering & Physics, 71:45–55, 2019.
  • [7] Daniel Canuto, Joe L Pantoja, Joyce Han, Erik P Dutson, and Jeff D Eldredge. An ensemble kalman filter approach to parameter estimation for patient-specific cardiovascular flow modeling. Theoretical and Computational Fluid Dynamics, pages 1–24, 2020.
  • [8] Caroline Cheng, Frank Helderman, Dennie Tempel, Dolf Segers, Beerend Hierck, Rob Poelmann, Arie van Tol, Dirk J Duncker, Danielle Robbers-Visser, Nicolette TC Ursem, et al. Large variations in absolute wall shear stress levels within one species and between species. Atherosclerosis, 195(2):225–235, 2007.
  • [9] Francisco J Contijoch, Michael Horowitz, Evan Masutani, Seth Kligerman, and Albert Hsiao. 4d flow vorticity visualization predicts regions of quantitative flow inconsistency for optimal blood flow measurement. Radiology: Cardiothoracic Imaging, 2(1):e190054, 2020.
  • [10] Luca Dede. Optimal flow control for navier–stokes equations: drag minimization. International Journal for Numerical Methods in Fluids, 55(4):347–366, 2007.
  • [11] Petter Dyverfeldt, Malenka Bissell, Alex J Barker, Ann F Bolger, Carl-Johan Carlhäll, Tino Ebbers, Christopher J Francios, Alex Frydrychowicz, Julia Geiger, Daniel Giese, et al. 4d flow cardiovascular magnetic resonance consensus statement. Journal of Cardiovascular Magnetic Resonance, 17(1):1–19, 2015.
  • [12] Marta D’Elia, Lucia Mirabella, Tiziano Passerini, Mauro Perego, Marina Piccinelli, Christian Vergara, and Alessandro Veneziani. Applications of variational data assimilation in computational hemodynamics. In Modeling of Physiological Flows, pages 363–394. Springer, 2012.
  • [13] Marta D’Elia and Alessandro Veneziani. Uncertainty quantification for data assimilation in a steady incompressible navier-stokes problem. ESAIM: Mathematical Modelling and Numerical Analysis, 47(4):1037–1057, 2013.
  • [14] Parastou Eslami, Justin Tran, Zexi Jin, Julia Karady, Romina Sotoodeh, Michael T Lu, Udo Hoffmann, and Alison Marsden. Effect of wall elasticity on hemodynamics and wall shear stress in patient-specific simulations in the coronary arteries. Journal of biomechanical engineering, 142(2), 2020.
  • [15] Mahdi Esmaily Moghadam, Yuri Bazilevs, Tain-Yen Hsia, Irene Vignon-Clementel, and Alison Marsden. A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations. Computational Mechanics, 48(3), 2011.
  • [16] Simon Wolfgang Funke, Magne Nordaas, Øyvind Evju, Martin Sandve Alnæs, and Kent Andre Mardal. Variational data assimilation for transient blood flow simulations: Cerebral aneurysms as an illustrative example. International Journal for Numerical Methods in Biomedical Engineering, 35(1):e3152, 2019.
  • [17] Xinyang Ge, Zhaofang Yin, Yuqi Fan, Yuri Vassilevski, and Fuyou Liang. A multi-scale model of the coronary circulation applied to investigate transmural myocardial flow. International Journal for Numerical Methods in Biomedical Engineering, 34(10):e3123, 2018.
  • [18] Frank JH Gijsen, Jolanda J Wentzel, Attila Thury, Frits Mastik, Johannes A Schaar, Johan CH Schuurbiers, Cornelis J Slager, Wim J van der Giessen, Pim J de Feyter, Anton FW van der Steen, et al. Strain distribution over plaques in human coronary arteries relates to shear stress. American Journal of Physiology-Heart and Circulatory Physiology, 295(4):H1608–H1614, 2008.
  • [19] Leopold Grinberg and George Em Karniadakis. Outflow boundary conditions for arterial networks with multiple outlets. Annals of biomedical engineering, 36(9):1496–1514, 2008.
  • [20] Harald C Groen, Lenette Simons, Quirijn JA van den Bouwhuijsen, E Marielle H Bosboom, Frank JH Gijsen, Alina G van der Giessen, Frans N van de Vosse, Albert Hofman, Antonius FW van der Steen, Jacqueline CM Witteman, et al. Mri-based quantification of outflow boundary conditions for computational fluid dynamics of stenosed human carotid arteries. Journal of Biomechanics, 43(12):2332–2338, 2010.
  • [21] Max D Gunzburger. Perspectives in flow control and optimization. SIAM, 2002.
  • [22] Michael Hinze, René Pinnau, Michael Ulbrich, and Stefan Ulbrich. Optimization with PDE constraints, volume 23. Springer Science & Business Media, 2008.
  • [23] Kazufumi Ito and Karl Kunisch. Lagrange multiplier approach to variational problems and applications. SIAM, 2008.
  • [24] Hyun Jin Kim, Irene E Vignon-Clementel, C Alberto Figueroa, John F LaDisa, Kenneth E Jansen, Jeffrey A Feinstein, and Charles A Taylor. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering, 37(11):2153–2169, 2009.
  • [25] Taha S Koltukluoğlu and Pablo J Blanco. Boundary control in computational haemodynamics. Journal of Fluid Mechanics, 847:329–364, 2018.
  • [26] Jeffrey C Lagarias, James A Reeds, Margaret H Wright, and Paul E Wright. Convergence properties of the nelder–mead simplex method in low dimensions. SIAM Journal on Optimization, 9(1):112–147, 1998.
  • [27] Rajnesh Lal, Bijan Mohammadi, and Franck Nicoud. Data assimilation for identification of cardiovascular network characteristics. International Journal for Numerical Methods in Biomedical Engineering, 33(5):e2824, 2017.
  • [28] WARREN K Laskey, H GRAHAM Parker, VICTOR A Ferrari, WILLIAM G Kussmaul, and ABRAHAM Noordergraaf. Estimation of total systemic arterial compliance in humans. Journal of Applied Physiology, 69(1):112–119, 1990.
  • [29] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.
  • [30] David M Martin, Eoin A Murphy, and Fergal J Boyle. Computational fluid dynamics analysis of balloon-expandable coronary stents: influence of stent and vessel deformation. Medical engineering & physics, 36(8):1047–1056, 2014.
  • [31] Massimiliano Mercuri. Tuning of boundary conditions parameters for hemodynamics simulation using patient data. PhD thesis, University of Sheffield, 2019.
  • [32] Umberto Morbiducci, Diego Gallo, Diana Massai, Filippo Consolo, Raffaele Ponzini, Luca Antiga, Cristina Bignardi, Marco A Deriu, and Alberto Redaelli. Outflow conditions for image-based hemodynamic models of the carotid bifurcation: implications for indicators of abnormal flow. Journal of Biomechanical Engineering, 132(9), 2010.
  • [33] Cecil D Murray. The physiological principle of minimum work: I. the vascular system and the cost of blood volume. Proceedings of the National Academy of Sciences of the United States of America, 12(3):207, 1926.
  • [34] Federico Negri, Andrea Manzoni, and Gianluigi Rozza. Reduced basis approximation of parametrized optimal flow control problems for the stokes equations. Computers & Mathematics with Applications, 69(4):319–336, 2015.
  • [35] Sanjay Pant, Benoit Fabrèges, J-F Gerbeau, and IE Vignon-Clementel. A methodological paradigm for patient-specific multi-scale cfd simulations: from clinical measurements to parameter estimates for individual analysis. International journal for numerical methods in biomedical engineering, 30(12):1614–1648, 2014.
  • [36] T. J. Pedley. The Fluid Mechanics of Large Blood Vessels. Cambridge University Press, 1980.
  • [37] Paris Perdikaris and George Em Karniadakis. Model inversion via multi-fidelity bayesian optimization: a new paradigm for parameter estimation in haemodynamics, and beyond. Journal of The Royal Society Interface, 13(118):20151107, 2016.
  • [38] S Pirola, Z Cheng, OA Jarral, DP O’Regan, JR Pepper, T Athanasiou, and XY Xu. On the choice of outlet boundary conditions for patient-specific analysis of aortic flow using computational fluid dynamics. Journal of biomechanics, 60:15–21, 2017.
  • [39] Alfio Quarteroni and Silvia Quarteroni. Numerical models for differential problems, volume 2. Springer, 2009.
  • [40] Alfio Quarteroni and Gianluigi Rozza. Optimal control and shape optimization of aorto-coronaric bypass anastomoses. Mathematical Models and Methods in Applied Sciences, 13(12):1801–1823, 2003.
  • [41] Robert S Reneman and Arnold PG Hoeks. Wall shear stress as measured in vivo: consequences for the design of the arterial system. Medical & Biological Engineering & Computing, 46(5):499–507, 2008.
  • [42] Rodrigo M Romarowski, Adrien Lefieux, Simone Morganti, Alessandro Veneziani, and Ferdinando Auricchio. Patient-specific cfd modelling in the thoracic aorta with pc-mri–based boundary conditions: A least-square three-element windkessel approach. International Journal for Numerical Methods in Biomedical Engineering, 34(11):e3134, 2018.
  • [43] Sethuraman Sankaran, Mahdi Esmaily Moghadam, Andrew M Kahn, Elaine E Tseng, Julius M Guccione, and Alison L Marsden. Patient-specific multiscale modeling of blood flow for coronary artery bypass graft surgery. Annals of Biomedical Engineering, 40(10):2228–2242, 2012.
  • [44] Meena Sankaranarayanan, Leok Poh Chua, Dhanjoo N Ghista, and Yong Seng Tan. Computational model of blood flow in the aorto-coronary bypass graft. Biomedical Engineering online, 4(1):14, 2005.
  • [45] Daniele E Schiavazzi, Alessia Baretta, Giancarlo Pennati, Tain-Yen Hsia, and Alison L Marsden. Patient-specific parameter estimation in single-ventricle lumped circulation models under uncertainty. International Journal for Numerical Methods in Biomedical Engineering, 33(3):e02799, 2017.
  • [46] Hang Si. Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software (TOMS), 41(2):1–36, 2015.
  • [47] J Soulis, G Giannoglou, M Dimitrakopoulou, V Papaioannou, S Logothetides, and D Mikhailidis. Influence of oscillating flow on ldl transport and wall shear stress in the normal aortic arch. The Open Cardiovascular Medicine Journal, 3:128, 2009.
  • [48] Ryan L Spilker and Charles A Taylor. Tuning multidomain hemodynamic simulations to match physiological measurements. Annals of Biomedical Engineering, 38(8):2635–2648, 2010.
  • [49] Zoran Stankovic, Bradley D Allen, Julio Garcia, Kelly B Jarvis, and Michael Markl. 4d flow imaging with mri. Cardiovascular Diagnosis and Therapy, 4(2):173, 2014.
  • [50] Martin Stoll and Andy Wathen. All-at-once solution of time-dependent stokes control. Journal of Computational Physics, 232(1):498–515, 2013.
  • [51] Charles A Taylor and CA Figueroa. Patient-specific modeling of cardiovascular mechanics. Annual review of biomedical engineering, 11:109–134, 2009.
  • [52] Charles A Taylor, Timothy A Fonte, and James K Min. Computational fluid dynamics applied to cardiac computed tomography for noninvasive quantification of fractional flow reserve: scientific basis. Journal of the American College of Cardiology, 61(22):2233–2241, 2013.
  • [53] Jorge Tiago, Alberto Gambaruto, and Adelia Sequeira. Patient-specific blood flow simulations: setting dirichlet boundary conditions for minimal error with respect to measured data. Mathematical modelling of natural phenomena, 9(6):98–116, 2014.
  • [54] Jorge Tiago, Telma Guerra, and Adélia Sequeira. A velocity tracking approach for the data assimilation problem in blood flow simulations. International Journal for Numerical Methods in Biomedical Engineering, 33(10):e2856, 2017.
  • [55] Bram Trachet, Joris Bols, Gianluca De Santis, Stefaan Vandenberghe, Bart Loeys, and Patrick Segers. The impact of simplified boundary conditions and aortic arch inclusion on cfd simulations in the mouse aorta: a comparison with mouse-specific reference data. Journal of Biomechanical Engineering, 133(12), 2011.
  • [56] Justin S Tran, Daniele E Schiavazzi, Abhay B Ramachandra, Andrew M Kahn, and Alison L Marsden. Automated tuning for parameter identification and uncertainty quantification in multi-scale coronary simulations. Computers & fluids, 142:128–138, 2017.
  • [57] Adam Updegrove, Nathan M Wilson, Jameson Merkow, Hongzhi Lan, Alison L Marsden, and Shawn C Shadden. Simvascular: an open source pipeline for cardiovascular simulation. Annals of Biomedical Engineering, 45(3):525–541, 2017.
  • [58] Alina G Van der Giessen, Harald C Groen, Pierre-Andre Doriot, Pim J De Feyter, Antonius FW Van der Steen, Frans N Van de Vosse, Jolanda J Wentzel, and Frank JH Gijsen. The influence of boundary conditions on wall shear stress distribution in patients specific coronary trees. Journal of Biomechanics, 44(6):1089–1095, 2011.
  • [59] Irene E Vignon-Clementel, C Alberto Figueroa, Kenneth E Jansen, and Charles A Taylor. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Computer Methods in Applied Mechanics and Engineering, 195(29-32):3776–3796, 2006.
  • [60] Peter E Vincent, AM Plata, AAE Hunt, PD Weinberg, and SJ Sherwin. Blood flow in the rabbit aortic arch and descending thoracic aorta. Journal of The Royal Society Interface, 8(65):1708–1719, 2011.
  • [61] Nico Westerhof, Jan-Willem Lankhaar, and Berend E Westerhof. The arterial windkessel. Medical & biological engineering & computing, 47(2):131–141, 2009.
  • [62] Zakia Zainib, Francesco Ballarin, Stephen Fremes, Piero Triverio, Laura Jiménez-Juan, and Gianluigi Rozza. Reduced order methods for parametric optimal flow control in coronary bypass grafts, towards patient-specific data assimilation. International Journal for Numerical Methods in Biomedical Engineering, page e3367, 2019.
  • [63] Mair Zamir, Paula Sinclair, and Thomas H Wonnacott. Relation between diameter and flow in major branches of the arch of the aorta. Journal of Biomechanics, 25(11):1303–1310, 1992.