1. Introduction and motivation
In recent years, large improvements have been made in developing patientspecific computational models to predict blood flow in patients affected by various cardiovascular diseases [51, 18]. The availability of data obtained from noninvasive 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 NavierStokes 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, NavierStokes 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 supraaortic branches: brachiocephalic artery (BCA), left common carotid artery (LCC), and left subclavian artery (LSUB).
The selection of appropriate boundary conditions is always a crucial step when dealing with patientspecific 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 patientspecific 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 (PCMRI), a technique for measuring blood velocity in specific sections of the vessel of interest, or with 4DFlow MRI, a timeresolved PCMRI 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 patientspecific 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 zerodimensional (lumped parameter) models[44, 43, 52, 17, 19], which, by prescribing a specific pressureflow relationship at each outlet, lead to more physiological pressure and flow waveforms, and to a better approximation of invivo data. Lumped models may consist of single resistances or more sophisticated zerodimensional models, such as the threeelement 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 patientspecific invivo 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 patientspecific data through a manual tuning process, until the userdesired accuracy is obtained [5, 31]. Despite being straightforward to implement, this approach is suboptimal, operatordependent, timeconsuming, and affected by nonrepeatability. To achieve robust and reliable patientspecific 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 nonlinear extensions, which performs the assimilation of measurements
on the fly and updates the unknown parameters accordingly [27, 7]. Arthurs et al. [4] proposed a reducedorder 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 patientspecific 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 hardtoquantify 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 steadystate blood flow simulations was presented. The methodology was validated against real 4Dflow 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 PCMRI data as outlet boundary conditions, tuning the parameters of a threeelement Windkessel model via a leastsquares approach. Zainib et al. [62] proposed a reduced order framework for the application of optimal control to coronary artery bypass grafts, where syntheticallygenerated measurements were used to tune Neumanntype 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 resistivetype outlet boundary conditions based on the solution of an optimal control problem, which assimilates a set of patientspecific measurements, e.g., through 4DFlow 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 steadystate NavierStokes 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 VignonClementel et al.[59]. To ensure a better match between simulation results and invivo data, the inlet flow waveform measured with 4DFlow MRI is imposed as a patientspecific 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 supraaortic branches, assimilating real 4DFlow 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 timeaveraged 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 patientspecific 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 patientspecific cardiovascular configurations from clinically acquired images, and the acquisition of patientspecific 4DFlow MRI data. The proposed framework is represented in Fig. 2.
2.1. Data acquisition and anatomical reconstruction
The proposed methodology has been tested in four clinical cases from a singlecenter 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 noninvasively. 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 supraaortic branches was acquired using 320detector row CT scanner (Aquilion One, Canon Medical Systems). From CT images, the vessels surface was reconstructed using the opensource 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 invivo using a 4Dflow 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 retrogating and adaptive navigator respiratory gating. Imaging parameters were as follows: velocity encoding=150 cm/s, field of view=200420 mm x 248368 mm, spatial resolution=1.93.5 x 2.03.2 x 1.83.5 mm, temporal resolution=39.947.2 ms, flip angle=8. After 4DFlow MRI, diastolic blood pressure and systolic blood pressure were measured with brachial cuffbased 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 NavierStokes 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 NavierStokes model. This analysis is reported in Section 3.
Considering the large diameter of the aorta and supraaortic 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 4DFlow MRI data. A noslip condition is imposed at the vessel walls , assumed to be rigid and nonpermeable, while at the outlets a resistive type boundary condition is imposed, following the coupled multidomain method proposed by VignonClementel 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 4DFlow 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.
To solve the optimal control problem, we adopt the adjointbased 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 socalled 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 socalled coupled optimality system, which we solve through a oneshot approach [21, 50], where the system is solved directly for all the unknown variables.
We adopt an optimisethendiscretise 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 finitedimensional solution spaces . In particular, we use TaylorHood elements for the velocitypressure pair, i.e. finite elements to define and finite elements for . Since the spaces associated to control and Lagrange multipliers are already finitedimensional, we set and . The discretised system is solved using the opensource libraries FEniCS [2, 29] and multiphenics [1], the latter being an opensource 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 patientspecific 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, threeelements 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 
3.2. Patientspecific data
Moving to the real cases with patientspecific 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 crosssectional 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 noninvasively on the patient, as reported in Section 2.1, divided by the mean aortic flow rate measured with 4DFlow 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 patientspecific measurements to compute flow splitting: instead, assuming that branches with larger crosssection have higher flow rates, it is solely based on the patient’s anatomy. This feature justifies its use in those studies where invivo 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 noninvasively and the outlet flow rates measured with 4DFlow 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 invivo. 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 derivativefree 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 3Dmodel of the underlying system, whereas here a 0D approximation is used.
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 4DFlow MRI data could be due to the finite resolution of 4Dflow 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 4DFlow 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 noninvasively after MRI, and the flow rates measured invivo with 4DFlow 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 nonphysical. 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.






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 % 

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 
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 4DFlow 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).
3.5. Unsteady NavierStokes simulations
An accurate representation of blood flow in the aorta is generally obtained through unsteady NavierStokes simulations, which provide a more realistic and accurate time evolution of blood flow. Assuming laminar regime, the incompressible NavierStokes equations take the form
(15) 
We used the resistance values reported in Table 4 as outlet boundary conditions of highfidelity unsteady NavierStokes 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 nonlinear, unsteady scenario. Second, it allows to analyse the impact that boundary conditions obtained with different estimation techniques have on wall shear stressrelated indicators, which are clinically relevant. For unsteady NavierStokes simulation, the type of outlet boundary condition which best replicates realistic flow conditions is the threeelement 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 timestep 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 NavierStokes 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 NavierStokes one, mostly due to the nonlinearity and timedependency introduced by the latter. Nevertheless, similar trends to those of the Stokes experiments reported in Table 3 can be observed when moving to NavierStokes 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 highfidelity NavierStokes simulations, thus supporting the approach of estimating BCs on a linearised Stokes model.
Additionally, Figure 24 reports a comparison of timedependent flow rates extracted from 4DFlow MRI and those obtained from timedependent simulations with boundary conditions estimated with optimal control. Even if the proposed estimation method assimilates only the average flow rates, the timedependent 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.
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 timeaveraged wall shear stress (TAWSS) and the oscillatory shear index (OSI), calculated from NavierStokes 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 ‘hotspots’, 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 operatorindependent technique for boundary condition estimation.
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 invivo 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 invivo measurements can be either pressure or flow waveforms. In particular, the proposed method does not specifically require data from 4DFlow 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 PCMRI modality. Moreover, the proposed formulation assimilates data in a leastsquares 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 steadystate, linearized NavierStokes 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 timedependent optimal control problem is at the moment computationally intractable. The necessity to use a steadystate formulation of the problem requires resorting to the Stokes model, as NavierStokes equations may not converge to a steadystate 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 timedependent simulation is still necessary. One main consequence of using a steadystate formulation is the possibility to estimate only the total value of the resistance. Estimating the outlet capacitance, in fact, would require solving a timedependent NavierStokes 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 3element 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 timedependent 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 invivo 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 resistancetype boundary conditions, while assimilating invivo pressure and flow rate measurements. The experiments conducted on four patient anatomies revealed the validity of the proposed optimal controlbased technique for assimilating 4DMRI 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 stressrelated 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 interoperator 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 resistivetype 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 AROMACFD 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, JeanYves 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 HomerVanniasinkam, Stavroula Balabani, and Vanessa DíazZuccarini. Patientspecific 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 patientspecific 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 RobbersVisser, 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, CarlJohan 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 navierstokes 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 patientspecific simulations in the coronary arteries. Journal of biomechanical engineering, 142(2), 2020.
 [15] Mahdi Esmaily Moghadam, Yuri Bazilevs, TainYen Hsia, Irene VignonClementel, 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 multiscale 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 PhysiologyHeart 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. Mribased 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 VignonClementel, 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 threedimensional 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, KentAndre 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 balloonexpandable 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 imagebased 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, JF Gerbeau, and IE VignonClementel. A methodological paradigm for patientspecific multiscale 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 multifidelity 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 patientspecific 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 aortocoronaric 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. Patientspecific cfd modelling in the thoracic aorta with pcmri–based boundary conditions: A leastsquare threeelement 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. Patientspecific 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 aortocoronary bypass graft. Biomedical Engineering online, 4(1):14, 2005.
 [45] Daniele E Schiavazzi, Alessia Baretta, Giancarlo Pennati, TainYen Hsia, and Alison L Marsden. Patientspecific parameter estimation in singleventricle lumped circulation models under uncertainty. International Journal for Numerical Methods in Biomedical Engineering, 33(3):e02799, 2017.
 [46] Hang Si. Tetgen, a delaunaybased 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. Allatonce solution of timedependent stokes control. Journal of Computational Physics, 232(1):498–515, 2013.
 [51] Charles A Taylor and CA Figueroa. Patientspecific 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. Patientspecific 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 mousespecific 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 multiscale 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, PierreAndre 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 VignonClementel, C Alberto Figueroa, Kenneth E Jansen, and Charles A Taylor. Outflow boundary conditions for threedimensional finite element modeling of blood flow and pressure in arteries. Computer Methods in Applied Mechanics and Engineering, 195(2932):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, JanWillem 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énezJuan, and Gianluigi Rozza. Reduced order methods for parametric optimal flow control in coronary bypass grafts, towards patientspecific 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.
Comments
There are no comments yet.