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

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 (WSS), 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.

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. , 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 Figure 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). 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). , , 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. 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. 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. A better solution resorts to the use of zero‐dimensional (lumped parameter) models, , , , , 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. 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.

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 standard approach for cardiovascular CFD simulations consists in a first estimation of outlet resistive BCs through Murray's law, 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. , 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. , 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. , Arthurs et al. 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. , , , Schiavazzi et al., 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. proposed an iterative approach based on adaptive Markov chain Monte Carlo sampling to quantify confidence in local hemodynamic indicators, for example, WSS 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. , , , 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. , , For what concerns the calibration of hard‐to‐quantify boundary conditions, Tiago et al., 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., 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. 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. 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., 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, WSS 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, for example, 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. 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 Figure 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 WSS and OSI 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.

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 Figure 2.

FIGURE 2

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

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 3 and 6 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. The reconstructed volume was discretised into tetrahedral elements using TetGen. After the CT scan, the blood velocity in the aorta and its branches was acquired in‐vivo using a 4D‐flow MRI sequence, using a 3 T 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 × 248–368 mm, spatial resolution = 1.9–3.5 × 2.0–3.2 × 1.8–3.5 mm3, temporal resolution = 39.9–47.2 ms, flip angle = 8°. After 4D‐Flow MRI, diastolic blood pressure P and systolic blood pressure P were measured with brachial cuff‐based method. The mean arterial pressure P was computed as

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 Figure 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 1 ≤ i ≤ i (i  = 4 in the case of Figure 1), respectively. The control variables are denoted as R . Thus, state equations can be written in strong form as where v is blood velocity, p the pressure, ν = .04 dynes/cm2s the dynamic viscosity, and n 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. Since the optimal control framework for partial differential Equations is typically derived starting from the weak formulation of the state equations, we introduce Hilbert spaces V (Ω) and P(Ω) for the velocity v and the pressure p, respectively. In particular, we choose and P(Ω) = L 2(Ω). Furthermore, we denote by the space associated to the controls . We note that V (Ω) and P(Ω) are function spaces (i.e., v(x) and p(x) are functions depending on the spatial coordinate x ∈Ω), while U is simply an Euclidean space (i.e., R are scalar numbers, and not functions).

Starting from the strong form in Equation (2), the weak formulation is derived as: given R ∈ U, find such that for every where w and q are the test functions associated to velocity and pressure, respectively.

We now reformulate the term in Equation (3) to obtain an equivalent weak formulation which is more suitable for the forthcoming finite element discretization. We introduce a set of Lagrange multipliers λ , 1 ≤ i ≤ i , defined as

We further denote by the vector collecting the Lagrange multipliers, and Z its associated space. Therefore, the equivalent weak formulation is: given R ∈ U, find v ∈ V (Ω), p ∈ P(Ω), λ ∈ Z such that for every , where η collects the test “functions” associated to the Lagrange multipliers λ. This is the final system representing how the control R 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 J with the following form

The first term in Equation (5) represents the normalized difference between the state pressure p and the patient's average pressure p measured at a specific cross section Γ ; as explained in Section 2.1, in this work p 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 Q 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:

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

To solve the optimal control problem, we adopt the adjoint‐based Lagrangian approach. , , , 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, z ∈ V 0(Ω), b ∈ P(Ω) and . The Lagrangian functional for this problem takes the form Given Equation (6), the optimality system can be obtained by imposing that the derivatives of with respect to (v, p, λ, R, z, b, t) must vanish, in short requiring . Taking the derivative of Equation (6) with respect to v (denoted by ) in the direction w we obtain while taking the derivative of Equation (6) with respect to p in the direction q we get

Similarly, the derivatives of Equation (6) with respect to λ , R , z, b, t are, respectively,

The presence in Equation (7) 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 Equation (7), and we add the following equations to the system

Equation (7) through Equation (15) form the so‐called coupled optimality system, which we solve through a one‐shot approach, , where the system is solved directly for all the unknown variables.

We adopt an optimize‐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 V (Ω), P (Ω), U , Z . In particular, we use Taylor‐Hood elements for the velocity–pressure pair, that is, finite elements to define V (Ω) and finite elements for P (Ω). Since the spaces associated to control and Lagrange multipliers are already finite‐dimensional, we set U  = U and Z  = Z. The discretised system is solved using the open‐source libraries FEniCS , and multiphenics, 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, a parallel sparse direct solver.

NUMERICAL RESULTS

Synthetic data

A mesh convergence analysis was first conducted on the patient‐specific meshes used for the experiments. In Figure 3 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 2.5 × 106 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.

FIGURE 3

Mesh convergence analysis

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 Q = 119.1 cm3/s, 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, R  = .09R , and a distal one, R  = .91R , as suggested by Kim et al. For the capacitance, a total value of .001 cm5/dyn was assumed, 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 Figure 4, which shows how the resistances chosen with optimal control allow to reconstruct the original flow waveform with a relative error of .09% for BCA, .09% for LCC, .1% for LSUB and .02% for DAo. The pressure waveform is recovered with a relative error of .005%.

TABLE 1

Results of experiment with synthetic data

	Resistance (dyn·s/cm5)				Flow rate (cm3/s)
	BCA	LCC	LSUB	DAo	BCA	LCC	LSUB	DAo
Ground truth estimated	7000	21,000	16,000	1700	19.02	6.30	8.13	80.02
	6937	20,846	15,991	1685	19.05	6.30	8.09	80.03

FIGURE 4

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

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, 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 where Q is the flow rate at the outlet, r is the radius at the outlet, and n for the aortic arch is conventionally set to 2. Multiplying r 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 where is the sum of the area of all the aortic outlets, while is the area of the outlet to which resistance R is associated. The total resistance R was computed as the mean pressure p measured non‐invasively on the patient, as reported in Section 2.1, divided by the mean aortic flow rate Q 0 measured with 4D‐Flow MRI, and then split among the outlet branches according to Equation (17). 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, , , , 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. 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. , ,

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 p computed from diastolic and systolic pressure measured non‐invasively and the outlet flow rates Q measured with 4D‐Flow MRI, outlet resistances can be computed as

This method is also common and it is based on the idea of performing the parameter estimation on 0D models 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 15% deviation, due to measurement uncertainty. For this reason, we developed a minimization problem which, using Ohm's law, estimates outlet resistances R 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 Equation (18). This is achieved by approximating the aorta with its equivalent 0D model, which is represented in Figure 5(A), and minimizing the cost function

FIGURE 5

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

The cost function was minimized using the Matlab function fminsearch, which employs a derivative‐free simplex search method. The cost functional reported in Equation (19) 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 (p ), and the second the outlet flow rates Q . However, optimal control relies on Stokes equations as a 3D‐model 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, 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. The amount of net flow is below 15% for all the cases under analysis, which is considered acceptable for 4D‐Flow MRI measurements according to Dyverfeldt et al. While in case 1 and case 4 the net flow is significant, for cases 2 and 3 it is practically negligible.

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

Case number	Inlet flow (cm3/s)	Outlet flow (cm3/s)	Net flow (cm3/s)	Mass conservation violation (%)
1	119.10	103.46	15.64	13%
2	107.00	107.18	−.18	.17%
3	125.63	125.60	.03	.02%
4	103.00	90.21	12.79	12.4%

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 (Equation (18)) are reported in the second row, while the third row contains the results obtained with the optimization based on Ohm's law as in Equation (19). 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 1% 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 min (wall clock CPU time), running on 18 Lenovo SD530 nodes, each with 40 Intel “Skylake” cores and 202 GB RAM.

TABLE 3

Comparison of pressure and flow rates for Stokes simulations with boundary conditions obtained with the different methods

Case number	Method	Pressure (mmHg)	Flow rate (cm3/s) (% error w.r.t. measurements)
			BCA	LCC	LSUB	DAo
1	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 (.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%)
2	Measurements	105	13.2	6.71	7.47	79.8
	Ohm's law	105 (0%)	13.1 (−.7%)	6.68 (−.4%)	7.44 (−.4%)	79.5 (−.4%)
	Opt. based on Ohm's law	105 (0%)	13.1 (−.7%)	6.70 (−.14%)	7.45 (−.26%)	79.6 (−.3%)
	Murray's law	106 (.9%)	17.0 (29%)	6.28 (6%)	9.30 (24%)	74.2 (−7%)
	Proposed	106 (.9%)	13.2 (0%)	6.71 (0%)	7.47 (0%)	79.5 (−.4%)
3	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 (.1%)
	Opt. based on Ohm's law	104 (.65%)	19.1 (.5%)	11.4 (.62%)	10.5 (0%)	84.9 (.1%)
	Murray's law	104 (.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 (.2%)
4	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%)

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 where u is a scalar control variable and v is the inlet velocity profile, with average value equal to the one measured with 4D‐Flow MRI. By choosing the best value for u , 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 Q at the inlet Γ :

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 (7)–(15). 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 min of wall clock CPU time, running on 24 Lenovo SD530 nodes, each with 40 Intel “Skylake” cores and 202 GB RAM).

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

Case number	Method	Pressure (mmHg)	Flow rate (cm3/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 (.5%)	8.53 (.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 (.8%)	4.33 (.23%)	6.96 (.6%)	73.0 (6%)

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

We used the resistance values R reported in Table 4 as outlet boundary conditions of high‐fidelity unsteady Navier–Stokes simulations performed using SimVascular. 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 analyze the impact that boundary conditions obtained with different estimation techniques have on WSS‐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 represented in Figure 5(B). Referring to Figure 5(B), each resistance R was split into a proximal one, R  = .09R , and a distal one, R  = .91R , as suggested by Kim et al. For the capacitance, a total value of .001 cm5/dyn was assumed, 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 ν = .04 dynes/cm2s, 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 .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 h (clock wall CPU time), running on 4 Lenovo SD530 nodes, each with 40 Intel “Skylake” cores and 202 GB RAM.

TABLE 4

Resistance values chosen by the different methods

Case number	Method	Resistance (dyn·s/cm5)
		BCA	LCC	LSUB	DAo
1	Ohm's law	8288	21,979	15,518	1800
	Opt. based on Ohm's law	8190	21,981	15,511	1679
	Murray's law	6837	21,242	17,591	1527
	Proposed	7941	21,609	15,153	1497
2	Ohm's law	10,690	21,003	18,847	1764
	Opt. based on Ohm's law	10,737	20,987	18,867	1767
	Murray's law	8306	22,451	15,161	1900
	Proposed	10,699	21,006	18,864	1773
3	Ohm's law	7248	12,142	13,094	1624
	Opt. based on Ohm's law	7255	12,145	13,090	1629
	Murray's law	7646	14,723	19,102	1513
	Proposed	7249	12,131	13,078	1621
4	Ohm's law	13,600	31,060	19,391	1943
	Opt. based on Ohm's law	13,500	31,069	19,399	1815
	Murray's law	8751	32,800	26,358	1708
	Proposed	13,166	30,500	18,903	1634

Figure 6 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 Figure 7. 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.

FIGURE 6

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

FIGURE 7

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

Additionally, Figure 8 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.

FIGURE 8

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 OSI, calculated from Navier–Stokes simulation results using the equations reported by Martin et al. For the sake of space, we report the analysis for case 3, but similar results were obtained for the other cases. Figure 9 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 Figure 10. For each point of the surface anatomy, the local relative difference was computed as:

FIGURE 9

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 10

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

For case 3, the difference in the TAWSS reaches a maximum relative difference of 24.8%, while the discrepancy in the OSI value reaches 55%. 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 (46%). 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.

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, 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 8, 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.

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 WSS‐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 modeling 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.