Physics-informed neural networks for blood flow inverse problems
Physics-informed neural networks (PINNs) have emerged as a powerful tool for
solving inverse problems, especially in cases where no complete information
about the system is known and scatter measurements are available. This is
especially useful in hemodynamics since the boundary information is often
difficult to model, and high-quality blood flow measurements are generally hard
to obtain. In this work, we use the PINNs methodology for estimating
reduced-order model parameters and the full velocity field from scatter 2D
noisy measurements in the ascending aorta. The results show stable and accurate
parameter estimations when using the method with simulated data, while the
velocity reconstruction shows dependence on the measurement quality and the
flow pattern complexity. The method allows for solving clinical-relevant
inverse problems in hemodynamics and complex coupled physical systems.
Jeremias Garay, Jocelyn Dunstan, Sergio Uribe, Francisco Sahli Costabal
# Physics-informed neural networks for blood flow inverse problems

###### Abstract

Physics-informed neural networks (PINNs) have emerged as a powerful tool for solving inverse problems, especially in cases where no complete information about the system is known and scatter measurements are available. This is especially useful in hemodynamics since the boundary information is often difficult to model, and high-quality blood flow measurements are generally hard to obtain. In this work, we use the PINNs methodology for estimating reduced-order model parameters and the full velocity field from scatter 2D noisy measurements in the ascending aorta. The results show stable and accurate parameter estimations when using the method with simulated data, while the velocity reconstruction shows dependence on the measurement quality and the flow pattern complexity. The method allows for solving clinical-relevant inverse problems in hemodynamics and complex coupled physical systems.

keywords: Physics-informed neural networks, hemodynamics, reduced-order modeling, blood flow, patient-specific model +
Footnote †: journal: Computer Methods in Applied Mechanics and Engineering

## 1 Introduction

Computational hemodynamics has been established as a relatively new research area in which blood flow is studied in many scenarios. Applications include the blood ejected from the heart through the aorta artery or blood flow through small vessels, such as the capillaries in the brain, or the flow through the coronaries arteries. Moreover, hemodynamics simulations have proven helpful for planning and assessing different cardiovascular pathologies, for instance, aortic dissection [1], stenosis [2], and aneurysms [3], which are known to be highly
affected by blood flow patterns and can generally be life-threatening and progressive. In a clinical context, the main reason for engaging in simulations in a complex cardiovascular scenario is to gain prediction capabilities, a feature that cannot be attained only based on images.

However, one of the main limitations of hemodynamical simulations and their clinical use is the strong sensitivity of the results to the modeling assumptions and parameters. Different imaging techniques have helped to fill this gap, helping to obtain patient-specific geometries and boundary condition information, which is essential for most physical models. Magnetic Resonance Imaging (MRI) [4] and Computational Tomography (CT) [5] are widely used and can be obtained relatively easily for small to large patient cohorts and healthy volunteers. Moreover, tissue and blood flow motion is usually measured by the so-called Phase-Contrast MRI (PC-MRI) technique [6], capable of encoding the tissue velocity into the phase of the emitted signal. From this method, we can obtain time-resolved vascular measurements, usually consisting of 20 to 30 snapshots within the cardiac cycle. Although the temporal and spatial resolution on different applications have become progressively better, some physical properties cannot be measured using purely non-invasive techniques, such as the stress at the boundary of the vessels, namely the wall shear stress, pressure drops across a stenotic region of the vessel, and other mechanical properties of the arteries.

Another limitation of hemodynamic simulations is the large computational cost and time required to make the models realistic. For this reason, several strategies for reducing model complexity have been proposed. For instance, the so-called _reduced order_ modeling approaches, which try to simplify some physical phenomena, making the simulation realistic while keeping a tractable computational time and needed resources [7; 8]. In hemodynamics, the range of the applications could vary in size, making the reduced strategy also change. Some examples on different applications are in: the coronary arteries [9; 10; 11], the aorta artery [12; 13], the pulmonary artery [14] and heart-valve problems [15; 16].

A popular choice for medium-to-large size arteries is using the three-element Windkessel model [17]. This model considers the vessels' compressibility and ability to store elastic energy during systole to release it in diastole. In general terms, a Windkessel model introduces a 0D differential equation at the vessel outlet, modeling the pressure and flow with an equivalent electric circuit, where the blood flow is represented by the current and the pressure as the voltage [18]. The model also introduces a resistance parameter, where higher resistance values can explain difficulties in the blood flow through that specific portion of the artery. Additionally, the model introduces a compliance \(C\), which considers the vessel's ability to store elastic energy.

The clinical value of estimating the Windkessel parameters on a specific data set relies on the physical information one can obtain from them, such as localized flow resistance effects or to know if a specific vessel presents a reduced elastic response. This data-driven information could be seen as a new cardiovascular biomarker for later clinical patient assessment due to its specificity and non-invasive nature. In this line, different strategies have been used to infer Windkessel parameters from partial measurements to make blood flow simulations patient-specific. For example, Arthurs et al. [19] applied the Reduced-Order Unscented Kalman Filter to the 3D time-dependent Navier-Stokes model with simplified fluid-solid
interaction effects on the wall to estimate the Windkessel parameters in the boundary conditions. Using variational data assimilation, Fevola et al. [20] estimated a single resistant model (i.e., \(C=0\)) on a stationary Stokes problem. Bertoglio et al. [21] used a Kalman Filter approach to estimate elastic and Windkessel parameters of the vessels from wall displacement measurements. Finally, Garay et al. [22] used a novel data assimilation term in their sequential inverse problem workflow to estimate three-element Windkessel parameters from highly aliased and noisy PC-MRI blood flow measurements.

However, all the works mentioned above strongly depend on the complete description of the physical model, sometimes leading to oversimplified solutions, especially in cases where the physics is complex and poorly understood. For that reason, in this work, we propose to use physics-informed neural networks (PINNs) [23; 24], which have shown great potential in solving inverse problems where the physical information is not complete, having the networks the ability to learn or to discover the missing parts from the data itself. We will use a feed-forward network architecture to represent the flow state of the coupled Navier-Stokes equations with the Windkessel model. Note that partial knowledge of the boundary conditions is very common in hemodynamics and usually leads to ill-posed problems which classical approaches failed to address [7]. Hence, PINNs allow us to estimate the time-dependent 3D velocity and pressure field from scatter-simulated medical images by solving the forward problem and the closest set of Windkessel parameters by solving the inverse problem simultaneously.

PINNs have emerged as a novel methodology for solving complex forward problems and ill-posed inverse problems in many research areas, such as in fluid mechanics [25; 26; 27], electrophysiology problems [28; 29], geosciences and wave-propagation problems [30; 31], non-linear solid mechanics [32; 33; 34] and also in cardiovascular biomechanics [35; 36; 37; 38]. See Ref [39] for an extensive review of all current applications. Although PINNs have been applied in many problems, and the Navier-Stokes equations have been focus of PINNs since its inception [23], to the best of the author's knowledge, this is the first time that PINNs have been applied in a coupled 3D fluid mechanics problem to estimate hemodynamic parameters from medical images.

The rest of this article is organized as follows. In Section 2, we present the mathematical model, the neural architecture to be used, and the reference velocity and simulated measurements used for the inverse problem. In Section 3, we present the results obtained for both studied regimes: the steady and transient flow cases. A discussion of the results and future lines of work are given in Section 4. Finally, the conclusions are presented in Section 5.

## 2 Methods

We start by introducing the model that describes the physics of blood flow in arteries.
### The mathematical model

Let \(\Omega\subset\mathbb{R}^{3}\) be the interior (lumen) of the thoracic aorta, represented in Figure 1, with its boundary \(\partial\Omega\) sub-divided as follows:

\[\partial\Omega=\Gamma_{in}\cup\Gamma_{wall}\cup\big{(}\cup_{k=1}\Gamma_{k}\big{)},\]

where \(\Gamma_{in}\) is the inlet boundary, \(\Gamma_{wall}\) the arterial wall, and \(\Gamma_{1},\ldots,\Gamma_{K}\) are the \(K\)-outlet boundaries. We consider that the blood flow in this domain is governed by the incompressible Navier-Stokes equations, with velocity \(\mathbf{u}(\mathbf{x},t)\) and pressure \(p(\mathbf{x},t)\):

\[\begin{cases}\rho\dfrac{\partial\mathbf{u}}{\partial t}+\rho\big{(}\mathbf{u} \cdot\nabla\big{)}\mathbf{u}-\mu\Delta\mathbf{u}+\nabla p=0\quad\text{in}\quad \Omega\times[0,T],\\ \nabla\cdot\mathbf{u}=0\quad\text{in}\quad\Omega\times[0,T],\\ \mathbf{u}=\mathbf{u}_{inlet}(t)\quad\text{on}\quad\Gamma_{in}\times[0,T],\\ \mathbf{u}=\mathbf{0}\quad\text{on}\quad\Gamma_{wall}\times[0,T],\\ \mu\dfrac{\partial\mathbf{u}}{\partial n}-p\mathbf{n}=-P_{k}(t)\mathbf{n} \quad\text{on}\quad\Gamma_{k}\times[0,T],\quad k=1,\ldots,K\end{cases} \tag{1}\]

where \(\rho\) is the blood density and \(\mu\) the dynamic viscosity of the fluid. Note that a Dirichlet inflow and a non-slip boundary conditions are adopted at \(\Gamma_{in}\) and \(\Gamma_{wall}\), respectively. Moreover, the pressure at the outlet \(P_{k}(t)\) is assumed given by the _three-element Windkessel_ model:

\[\begin{cases}P_{k}=R_{p,k}\ Q_{k}+\pi_{k},\\ Q_{k}=\int_{\Gamma_{k}}\mathbf{u}\cdot\mathbf{n},\\ C_{k}\dfrac{d\pi_{k}}{dt}+\dfrac{\pi_{k}}{R_{d,k}}=Q_{k}.\end{cases} \tag{2}\]

In the Windkessel model, \(R_{p,k}\) and \(R_{d,k}\) represent the resistance of the vasculature proximal and distal to \(\Gamma_{k}\), respectively, while \(C_{k}\) is the compliance of the distal vessels. The exterior normal vector of the outlet is represented by \(\mathbf{n}\).

In the special case of a stationary flow, all time derivatives vanish from both models, resulting in the lower dimensional system an effective _one-element_ Windkessel boundary condition with no internal pressure variable \(\pi\). The equations of the Windkessel model are reduced to:

\[\begin{cases}P_{k}=(R_{p,k}+R_{d,k})\ Q_{k}=R_{t,k}\ Q_{k},\\ Q_{k}=\int_{\Gamma_{k}}\mathbf{u}\cdot\mathbf{n}.\end{cases} \tag{3}\]

with the distal and proximal resistances absorbed in a total resistance parameter \(R_{t}\).
### Physics-Informed Neural Networks (PINNs)

We will use PINNs to solve the inverse problem of inferring the Windkessel model parameters \(R_{p,k},R_{d,k},C_{k}\quad\forall k=1,...,K\), and the velocity and pressure fields in the aorta from partial velocity measurements \(\mathbf{u}_{meas}\), obtained from PC-MRI and an average pressure measurement in the outlets \(\bar{p}_{meas}\). We represent the solution of the coupled system in Equation 1 and 2 by a neural network consisting of \(n\) feed-forward layers, as shown in Figure 2:

\[\left(\mathbf{u}(\mathbf{x},t),\ p(\mathbf{x},t)\right)=\mathcal{N}\mathcal{N} (\mathbf{x},t;\mathbb{W},\mathbf{b}) \tag{4}\]

Where the neural network \(\mathcal{N}\mathcal{N}\) outputs the fields \(\mathbf{u}\) and \(p\), and takes as input the coordinate points in space and time, parametrized by weights (\(\mathbb{W}\)) and biases (\(\mathbf{b}\)). The values for \(\mathbb{W}\) and \(\mathbf{b}\) are obtained after training the network. To achieve this, we define a loss function \(\mathcal{L}_{tot}\) that encourages learning the measurements, the normalized physics equations (1)a, (1)b and, (2), and the non-slip boundary condition in (1)d:

\[\mathcal{L}_{tot}=\mathcal{L}_{NS}+\mathcal{L}_{WK}+\mathcal{L}_{BC}+ \mathcal{L}_{u,data}+\mathcal{L}_{p,data}+\mathcal{L}_{gradp}. \tag{5}\]

where the individual loss components are defined as:

\[\mathcal{L}_{NS} = \lambda_{NS}\ \left|\left|\frac{\partial\mathbf{u}}{\partial t}+ \left(\mathbf{u}\cdot\nabla\right)\mathbf{u}-\frac{1}{Re}\Delta\mathbf{u}+ \nabla p\right|\right|_{\Omega}+\lambda_{NS}\ \left|\left|\nabla\cdot\mathbf{u}\right|\right|_{\Omega}, \tag{6}\] \[\mathcal{L}_{WK} = \lambda_{WK}\ \sum_{k=1}^{K}\left|\left|p_{k}-R_{p,k}Q_{k}-\pi_{k} \right|\right|_{\Gamma_{k}}+\left|\left|C_{k}\frac{d\pi_{k}}{dt}+\frac{\pi_{k} }{R_{d,k}}-Q_{k}\right|\right|_{\Gamma_{k}},\] (7) \[\mathcal{L}_{BC} = \lambda_{BC}\ ||\ \mathbf{u}\ ||_{\Gamma_{wall}},\] (8) \[\mathcal{L}_{u,data} = \lambda_{data}\ ||\ \mathbf{u}-\mathbf{u}_{meas}\ ||_{\Omega_{meas}},\] (9) \[\mathcal{L}_{p,data} = \lambda_{data}\ ||\ \bar{p}-\bar{p}_{meas}\ ||_{\Omega_{meas}},\] (10) \[\mathcal{L}_{p,gradp} = \lambda_{gradp}\sum_{k=1}^{K}\left|\left|\ \nabla p\ \right|\right|_{\Gamma_{k}}. \tag{11}\]

Figure 1: Aortic domain \(\Omega\). The Windkessel outlets are defined at \(\Gamma_{1}\) to \(\Gamma_{5}\) consisting of 2 resistors (\(R_{d}\), \(R_{p}\)) and 1 capacitor \(C\) each.
The last added term (11) is included to promote a null pressure gradient at the Windkessel outlets since the coupled equations between the 3D and 0D systems assumed the pressure to be constant at the interface. In all cases, \(||\cdot||\) represents the \(L_{2}\) norm of the quantities. We make all equations dimensionless for the loss function to obtain normalized predictions of velocity and pressure, which improves the training procedure. We use the following normalization of the physical quantities:

\[\mathbf{x} \longrightarrow\mathbf{x}/L\] \[\mathbf{u} \longrightarrow\mathbf{u}/U\] \[t \longrightarrow t/(L/U)\] \[p \longrightarrow p/(\rho U^{2})\]

where \(U\) and \(L\) are a characteristic velocity and length scale of the system, respectively. Consequently, the Reynolds number \(Re\), defined as \(Re=\frac{\rho UL}{\mu}\), appears in the Laplacian term of Equation (6). The Windkessel system is also dimensionally reduced as:

\[\pi \longrightarrow\pi/(\rho U^{2})\] \[Q \longrightarrow Q/(UL^{2})\] \[R_{d,p} \longrightarrow R_{d,p}/(\rho U/L^{2})\] \[C \longrightarrow C/(L^{3}/\rho U^{2})\]

To fix the blood pressure level of the coupled system, we also give to the inverse problem the average pressure curve (obtained from the reference solution) defined as:

\[\bar{p}=\sum_{k=0}^{K}\frac{1}{area(\Gamma_{k})}\int_{\Gamma_{k}}pdS. \tag{12}\]

Furthermore, we interpolate the latter curve in time to match the MRI measurements timestep, so both measurements set (velocity and mean pressure) are consistent in time. Note that for the steady problem, this curve consists of a single value representing the systolic pressure of the artery. In a clinical scenario, this mean pressure curve could be approximated as the pressure measured using a sphygmomanometer or an oscillometric blood pressure monitor, generally placed at the patient's upper arm or wrist [40; 41; 42].

Concerning the collocation points in which the term \(\mathcal{L}_{NS}\) is enforced, we use \(403,223\) points corresponding to a tetrahedral generated mesh for later usage. For each point, we generate a random time within the cardiac cycle as input to the network in the transient problem. Also, training by mini-batch is performed using only around \(1\%\) of these points at every iteration.

Furthermore, the coupling between the models requires the computation of a boundary integral (the velocity flow \(Q_{k}\) at \(\Gamma_{k}\)). For that reason, all the points at the outlets have to be evaluated at the same time. To optimize the computational resources, we evaluate the velocity of every outlet point only at the measurement times (i.e., \(22\) time intervals within the cardiac cycle) and at the end of every epoch. With this set of points, the flow is estimated
using a finite element quadrature and interpolated in time using a cubic interpolator for later evaluation at any time within the cardiac cycle.

The transient problem presents its unique challenges compared to the steady problem. For this reason, we take into account some additional considerations that were not necessary for the steady problem:

* We use a vector potential representation of the velocity \(\mathbf{u}\) to hard-impose mass conservation of the solution. Consequently, the velocity is assumed to be originated from the curl of a vector potential \(\Phi(\mathbf{x},t)\)[43] as: \[\mathbf{u}(\mathbf{x},t)=\nabla\times\Phi(\mathbf{x},t).\] (13) With this definition, the velocity field \(vecu\) is always divergence free. We find that this change of variable improves the quality of the results and convergence of the method at the expense increasing the computational time.
* We tune each weight of the total loss function automatically following the algorithm presented in [44], using the \(L_{1}\) norm of the gradient of the physics loss. We perform the update every 10 epochs. The new values of the weights are exponentially averaged with their respective history following the rule: \[\lambda_{new,i}=(1-\alpha)\cdot\lambda_{old,i}+\alpha\cdot\frac{||\nabla \mathcal{L}_{phys}||_{1}}{\overline{\nabla\mathcal{L}_{i}}}\qquad\forall i=1,2,...,N_{L}\]

Figure 2: Feed-forward neural network architecture used for the transient problem. The network represents a vector potential \(\Phi\) and the pressure field. After physics-informed training, the output of the network is used for reconstruct the 3D velocity field and the Windkessel parameters of all the domain outlets.
being \(N_{L}\) the number of loss terms. In all cases, we set \(\alpha\) to \(0.1\).
* We divide the training into two stages. First, we randomly initialize the Windkessel model parameters \(R_{p,k},R_{d,k},C_{k}\) and fix them. Then, we train only the network parameters \(\mathbb{W},\mathbf{b}\). In a second stage, we allow the optimizer to also change the Windkensell model parameters. We found that splitting the training in this way leads to better results since the network had time first to adapt the representation of physical velocities and pressures and then to estimate the Windkessel parameters.
* We compute the compliance parameter \(C\) from the distal resistance \(R_{d}\) and an estimation of the decay time \(R_{d}\cdot C\). This time was obtained at the beginning of the training from the mean pressure curve used as an extra measurement. Furthermore, the decay time was computed by fitting an exponential curve from the valve closure time and the end of the diastole, an interval at which the Windkessel model assumes a decay behavior modulated by \(\sim\exp(-t/R_{d}C)\). The decay time was assumed constant for all outlets at the aorta.

Finally, we obtained the weights and biases by minimizing the total loss function as:

\[\mathbb{W},\mathbf{b}=\arg\min_{W,b}\ \mathcal{L}_{tot}(W,b). \tag{14}\]

Figure 2 shows a schematic of the final architecture used for the transient problem. The minimization was performed using the ADAM optimizer [45] with the _swish_ activation function [46]. We perform the estimation of the velocity and pressure fields and the estimation of the Windkessel parameters simultaneously. The Winkessel parameters start from a random seed between a predefined confidence interval between the half and the double of the respective reference value and then evolve with the optimization. Five independent realizations are performed in which the network weights are randomly initialized in each one of these. All cases are summarized in Table 1.

The training of all the experiments were done in an NVIDIA RTX-A6000. The codes were written in Python using the PyTorch library [47] and are available on GitHub at [https://github.com/yeyemedicen/PINNs-WK-MRI](https://github.com/yeyemedicen/PINNs-WK-MRI). All used hyperparameters are reported for both stationary and transient problems in Table 2.

The normalized velocity and length values were tuned by hand, keeping in mind that the normalized quantities entering the PINNs workflow have to be ideally not larger than the unity. In the case of the transient problem, the value of the normalized velocity had to be significantly reduced because of the lower velocities at the end of the cardiac cycle (diastole).

\begin{table}
\begin{tabular}{c c c} \hline \hline  & Steady Problem & Transient Problem \\ \hline Estimation & \(\mathbf{u}\), \(p\), \(\left\{R_{tot}\right\}_{k=1}^{5}\) & \(\mathbf{u}\), \(p\), \(\left\{R_{p},R_{d},C_{d}\right\}_{k=1}^{5}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Different quantities to estimate in the steady and transient problems.
### The Reference Solution

To create the reference solution, equations (1) and (2) are solved using the finite element method (FEM), discretizing in space using stabilized \(\mathbb{P}1/\mathbb{P}1\) elements, and in time using a backward Euler scheme. The system was solved using a non-incremental fractional step scheme detailed in Ref [22]. The initial conditions were set to:

\[\begin{cases}\mathbf{u}=\mathbf{0},\\ \pi_{k}=78.8\ \mathrm{mmHg}\quad\mathrm{for}\quad\mathrm{k}=1,\dots,\mathrm{K}. \end{cases}\]

The values of \(\pi_{k}\) correspond to approximately the periodic state of the 3D-0D system taken from Ref [48]. Finally, we assume a velocity profile at \(\Gamma_{\mathit{inlet}}\) as:

\[\mathbf{u}(\mathbf{x},t)_{\mathit{inlet}}=\mathbf{u}_{\mathit{stokes}}( \mathbf{x})f(t)\ \mathbf{n}, \tag{15}\]

where \(\mathbf{u}_{\mathit{stokes}}(\mathbf{x})\) is the solution of a Stokes problem to ensure a parabolic-shape solution already adapted to the domain as is shown in Figure 3(a). The function \(f(t)\) represents the time-dependency of the inflow velocity and is taken in such a way that the total flow through \(\Gamma_{in}\) follows the curve in Figure 3(b), obtained from the Vascular Model Repository database [48]. The whole boundary condition is then added into the variational form as a penalization term \(A_{\mathit{inlet}}\) defined as:

\[A_{\mathit{inlet}}=\gamma\int_{\Gamma_{\mathit{inlet}}}(\mathbf{u}-\mathbf{u }_{\mathit{inlet}})\cdot\mathbf{v} \tag{16}\]

\begin{table}
\begin{tabular}{c c c} \hline \hline Hyperparameter & Steady Problem & Transient Problem \\ \hline Batchsize & \(1580\) & \(3000\) \\ Learning Rate & \(10^{-3}\) & \(10^{-3}\) \\ Learning Rate of Parameters & \(10^{-2}\) & \(10^{-2}\) \\ Learning Rate Scheduler & Yes & Yes \\ Hidden layers & \(7\) & \(7\) \\ No. neurons per layer & \(220\) & \(220\) \\ Hidden layers (\(\pi\)) & \(-\) & \(6\) \\ No. neurons per layer (\(\pi\)) & \(-\) & \(10\) \\ Activation Function & swish & swish \\ Epochs & \(1250\) & \(120+1250\) \\ Potential Vector Representation & No & Yes \\ Normalized Length (\(L\)) & \(1\) cm & \(0.5\) cm \\ Normalized Velocity (\(U\)) & \(5000\) cm/s & \(120\) cm/s \\ Loss weights & Manual & Automatic \\ \hline \(\lambda_{\mathit{phys}}\) & \(1.5\) & \(-\) \\ \(\lambda_{\mathit{data},u}\) & \(1.0\) & \(-\) \\ \(\lambda_{\mathit{data},p}\) & \(1.0\) & \(-\) \\ \(\lambda_{\mathit{BC}}\) & \(6.0\) & \(-\) \\ \(\lambda_{\mathit{winldk}}\) & \(1.0\) & \(-\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Hyperparameters used for training the PINNs
with the parameter \(\gamma\) fixed in \(10^{5}\) gr/(cm\({}^{2}\cdot\) s).

The Windkessel parameters were also obtained from the aforementioned repository database. For the numerical values of these constants, see Table 3.

The computational domain corresponds to a thoracic aorta with a total stream-wise length of 10.8 cm and meshed with 558,112 tetrahedrons and 115,779 vertices. We assume blood has a density of 1.06 gr/cm\({}^{3}\) and a constant dynamic viscosity of 0.035 P. The time step was set to \(\tau=0.001\) s with a total run time of 0.66 s. We also added backflow stabilization in every system outlet [49].

Figure 4 shows the resulting velocity field at peak systole next to each outlet's flow rates and pressure curves.

### PC-MRI measurement generation

From the reference solution, we simulate a 2D phase-contrast magnetic resonance acquisition located at three axial planes covering the whole domain, as shown in Figure 5(a). Thus, the velocity field is first interpolated into a voxel-like slice mesh with a resolution of \(1.0\times 1.0\) mm\({}^{2}\) and downsampled in time from \(1\)\(ms\) to \(30\)\(ms\), following the typical timestep of a clinical acquisition. After that, two different complex magnetization vectors are defined as:

\[M_{meas}^{u} = M_{0}\ \exp(i\phi_{0}+i\pi u_{MRI}/venc), \tag{17}\] \[M_{meas}^{0} = M_{0}\ \exp(i\phi_{0}), \tag{18}\]

Figure 3: (a) Stokes-profile velocity at the inlet of the aorta. (b) Velocity flow at the inlet taken from the SimVascular dataset.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Parameter & \(\Gamma_{1}\) & \(\Gamma_{2}\) & \(\Gamma_{3}\) & \(\Gamma_{4}\) & \(\Gamma_{5}\) \\ \hline \(R_{p}\) (dyn\(\cdot\) s\(\cdot\) cm\({}^{-5}\)) & 713 & 713 & 602 & 689 & 98 \\ \hline \(R_{d}\) (dyn\(\cdot\) s\(\cdot\) cm\({}^{-5}\)) & 12023 & 12023 & 10143 & 11609 & 1650 \\ \hline \(C_{d}\) (dyn\({}^{-1}\cdot\) cm\({}^{5}\)) & 8.256e-5 & 8.256e-5 & 9.785e-5 & 8.55e-5 & 6.015e-4 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Three-element Windkessel parameters for every outlet.
where \(M_{meas}^{u}\) is the velocity-encoded magnetization and \(M_{meas}^{0}\) is an extra measurement usually done to capture the reference phase \(\phi_{0}\). The interpolated velocity is represented by \(u_{MRI}\) while \(M_{0}\) is assumed constant within the vessel's lumen. Moreover, Gaussian noise is added to the magnetization components, producing a fixed signal-to-noise ratio of 18 dB. Consequently, a non-Gaussian noise distribution is induced in the reconstructed velocity, as is generally observed in PC-MRI measurements [50]. The final velocity field is then reconstructed using the following equation:

\[\mathbf{u}_{meas}=venc\ \frac{\angle\big{(}M_{meas}^{u}/M_{meas}^{0}\big{)}}{ \pi}, \tag{19}\]

with the symbol \(\angle\) representing the angle of the complex quantity measured from \(-\pi\) to \(\pi\). The velocity encoding parameter (\(venc\)) is chosen as the 120% of the maximum velocity to ensure aliased-free measurements. Figure 5(b) shows the resulting slice measurements for the stationary problem. We perform a similar procedure for the transient flow solution, resulting in 22 phases along the cardiac cycle.

## 3 Results

### Steady problem

In this section, we present the results obtained for the stationary problem. For each experiment, statistics are presented of 5 independent experiment realizations.

Figure 6 shows the evolution of the estimated total resistance by every training epoch and Windkessel outlet. We normalize all values to ease the visualization of the parameter

Figure 4: Reference simulation results. The blood flow velocity at peak systole and its flow rates and average pressure curves are depicted in (a) and (b), respectively.
evolution. It can be seen that all parameters stabilize rapidly to their reference value, yet oscillations are present; they reduce at the last stages of the training. Also, not all parameters show the same convergence rate since, for instance, the outlet number 3 takes around one-third of the total training epochs to converge to its reference value, while the outlet number 5 needs to approximately double the iterations to have a same-accurate result.

Table 4 shows the final values of the estimation, the mean estimated flow, and pressure at every outlet. These quantities were computed from the mean of the velocities and pressures by every experiment realization. We have found that the PINNs represent the solution at

\begin{table}
\begin{tabular}{l c c c c c c c c c} \hline \hline  & \multicolumn{8}{c}{Steady Problem} \\ \cline{2-10}  & \multicolumn{2}{c}{\(Q\) (mL/s)} & \multicolumn{4}{c}{\(P\) (mmHg)} & \multicolumn{4}{c}{\(R_{tot}\) (\(10^{3}\) dyn\(\cdot\) s\(\cdot\) cm\({}^{-5}\))} \\ \cline{2-10} Bnd. & Ref. & Mean & Std & Ref & Mean & Std & Ref & Mean & Std \\ \hline \(\Gamma_{in}\) & -298.88 & -297.16 & 2.07 & 248.39 & 248.92 & 0.63 & \(-\) & \(-\) & \(-\) \\ \(\Gamma_{1}\) & 25.82 & 25.58 & 0.10 & 247.96 & 246.93 & 1.63 & 12.74 & 12.86 & 0.14 \\ \(\Gamma_{2}\) & 25.69 & 25.64 & 0.13 & 247.90 & 247.34 & 0.80 & 12.74 & 12.89 & 0.06 \\ \(\Gamma_{3}\) & 30.56 & 30.71 & 0.27 & 248.20 & 247.98 & 0.66 & 10.75 & 10.78 & 0.08 \\ \(\Gamma_{4}\) & 26.65 & 26.08 & 0.12 & 248.22 & 249.08 & 0.66 & 12.30 & 12.71 & 0.10 \\ \(\Gamma_{5}\) & 188.85 & 185.97 & 0.84 & 248.15 & 248.30 & 0.66 & 1.74 & 1.78 & 0.04 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Flow and mean pressure for each outlet of the aorta. The mean pressures were computed as \(P_{i}=\frac{1}{Area_{i}}\int_{\Gamma_{i}}pds\), where \(i\) is the corresponding outlet number. The mean and standard deviation values reported are computed from the results of every seed.

Figure 5: Synthetic PC-MRI measurements obtained from the simulation. The axial planes taken from the aorta are shown in (a), while the resulting PC-MRI velocity measurements are shown in from frontal view in (b).
the inlet remarkably well, even though no information was given at that boundary. Also, the estimated parameters match their corresponding reference values well despite the lack of volumetric data and the artifacts introduced to the slice measurements due to the added noise and spatiotemporal interpolation. Finally, Figure 7 shows the velocity streamlines computed from the inlet nodes of the mesh of both reference and mean estimation velocities. This figure shows that the PINNs' solution is similar to the reference velocity while in the vessel's lumen and moderately fails for the nodes closer to the walls.

### Transient problem

Figure 8 shows the parameter estimation evolution. As before, we normalize the parameters by their reference values and split them in the vertical axis to ease the visualization. From the plots, we can see that proximal resistances (\(R_{p}\)) show significantly more oscillation during the first epochs of the training than the distal resistances parameters (\(R_{d}\)), though both show relatively fast convergence to their respective reference values within the half of the training. However, some outlets (\(\Gamma_{3}\) and \(\Gamma_{4}\)) present a growing standard deviation in the last epochs of the training for the estimated resistance parameters. Table 5 reports the final values with their respective statistics.

Figure 9 shows the obtained flows and mean pressure curves. From these curves, a general flow underestimation is shown in the upper outlets of the aorta within the systole, while the inlet and descending aorta (\(\Gamma_{in}\) and \(\Gamma_{5}\)) were remarkably well captured in the whole cardiac cycle. Conversely, mean pressure curves computed as:

\[P_{i}=\frac{1}{area(\Gamma_{i})}\int_{\Gamma_{i}}p\ dS,\]

w

Figure 6: Evolution of Parameter estimation statistics for each total resistances of the aorta in steady flow regime. All values were divided by their reference value and split vertically for better visualization. Dashed black lines represent the exact solution, and continuous color lines represent the estimation mean. The colored area surrounding the curves is determined by the minimum and maximum value found along the multiple realizations of the experiment
tend to be overestimated in the upper aorta and underestimated in the downward portion of the domain. Moreover, an oscillating behavior was found in the flows and pressure in the outlets \(\Gamma_{3}\) and \(\Gamma_{4}\), outlets where the standard deviation of the estimated parameters is the most. A FEM simulation was also performed after the training, using the estimated parameters and inlet velocity profile found in the experiment, which present the median error among all realizations. This solution could be seen as a post-process step of the PINNs' estimation and generally presents reduced oscillations compared with the raw estimated velocity and pressure fields. The obtained curves of the post-process FEM simulation are depicted in colored dashed lines.

Finally, streamlines computed from the nodes of \(\Gamma_{in}\) were calculated at two instants of the cardiac cycle and shown in Figure 10. This was done for the reference velocity, the PINN's estimation, and the FEM post-process simulated velocity. Generally, the neural network

\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline  & \multicolumn{8}{c}{Transient Problem} \\ \cline{2-9}  & \(R_{p}\) (\(10^{3}\) dyn \(\cdot\) s \(\cdot\) cm\({}^{-5}\)) & \multicolumn{2}{c}{\(R_{d}\) (\(10^{3}\) dyn \(\cdot\) s \(\cdot\) cm\({}^{-5}\))} & \multicolumn{2}{c}{\(C_{d}\) (\(10^{-5}\) dyn\({}^{-1}\cdot\) cm\({}^{5}\))} \\ \cline{2-9} Bnd. & Ref. & Mean & Std & Ref & Mean & Std & Ref & Mean & Std \\ \hline \(\Gamma_{1}\) & 0.71 & 1.02 & 0.10 & 12.02 & 15.76 & 1.87 & 8.26 & 6.42 & 0.70 \\ \(\Gamma_{2}\) & 0.71 & 1.10 & 0.12 & 12.02 & 16.75 & 1.58 & 8.26 & 6.01 & 0.52 \\ \(\Gamma_{3}\) & 0.60 & 0.91 & 0.15 & 10.14 & 14.51 & 1.88 & 9.78 & 6.99 & 0.82 \\ \(\Gamma_{4}\) & 0.69 & 0.84 & 0.14 & 11.61 & 14.46 & 1.94 & 8.55 & 7.01 & 0.80 \\ \(\Gamma_{5}\) & 0.10 & 0.12 & 0.02 & 1.65 & 2.02 & 0.22 & 60.15 & 50.18 & 5.51 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Parameter estimation obtained on the transient-flow regime.

Figure 7: Velocity streamlines when the steady-state flow regime is imposed for (a) the reference solution and (b) the mean estimated velocity represented by the PINNs
can capture the main flow features at systole but not so well during diastole, where the velocities are lower and the measurements' velocity-to-noise ratio is worse. When comparing the reference velocity against the FEM post-process solution, no significant gain in quality can be seen while in systole. However, during diastole, the FEM solution captured the main flow features that the PINNs' estimation failed to represent, especially the vortices at the ascending aorta and the near-wall velocity in the entire domain.

## 4 Discussion

In this we work, we develop a novel method to estimate hemodynamical parameters from sparse measurements using PINNs. We note that the original problem was divided into two different flow regimes: stationary and transient. We made this split to study the quality of the results when the problem grows in complexity.

The stationary problem could be seen as a inexpensive version of the transient problem, and in clinical hemodynamics, is usually studied when the transient problem's computational time exceeds the allowed limit imposed by the clinical application itself [51; 52]. In terms of the quality of the estimations, it could be seen from the results in Table 4 that although all parameters were overestimated, the obtained values present a mean error of less than 2% when comparing with their reference. On the other hand, when looking into the reconstructed velocity streamlines in Figure 7, both reconstructed and estimated fields are very similar, with some mild differences at the vessel wall. One of the reasons for that difference is because FEM methods exactly imposes the boundary condition on the walls, being the velocity there set exactly to zero, as it happened with the reference solution.

Figure 8: Evolution of the parameter estimation statistics for each Windkessel outlet of the aorta in transient flow regime. All values were divided by their reference value and split vertically for better visualization. Dashed black lines represent the exact solution, and continuous color lines represent the estimation mean. The colored area surrounding the curves is determined by the minimum and maximum values found along the multiple realizations of the experiment.
However, in the case of the PINNs method presented here, the non-slip condition enters the optimization as a loss term and, consequently, is not completely satisfied. This produces some of the velocity lines to go out of the vessel through the walls. This effect is also present in the transient flow regime.

It is particularly interesting that while the velocity tends to be misrepresented near the vessel wall, the outflow values are not. We hypothesized that since the coupling between the Navier-Stokes and Windkessel model is in terms of pressures and flows [17], the extra flow condition in the total loss in Equation (7) acts as a reinforcement for the training. Consequently, the outlet velocity tends to rise in the lumen, compensating for the low velocity near the wall. This effect is increased during the transient diastole, as some of the upper aortic outlets do not perceive any velocity, as shown in Figure 9(b,e). One strategy to mitigate this effect could be introducing a hard-constrain for the non-slip condition into the optimization, as done for a 2D Navier-Stokes system in Ref [53]. Another strategy could be relaxing the wall definition by introducing a _slip/transpiration_ condition as it was done

Figure 9: Flows (red) and mean pressure (sky blue) curves at every aortic outlet over time. Continuous lines represent mean estimation, while the area plot is defined by the minimum and maximum value encountered in every experiment realization. Black dashed lines correspond to the respective reference value, while colored dashed lines correspond to the curves obtained from the post-process FEM simulation.
in Nolte et al. [54] when studying the impact of segmentation's errors in MRI blood flow measurements.

Another aspect to discuss more in detail is the use of the potential vector representation for the flow velocity in Equation (13). This strategy could be seen as a hard constraint, applied only to the transient problem, to exactly satisfy the mass-conservation term of the Navier-Stokes model, as it was also done in Ref [55] for turbulent flows. We found that this change of variables improved the overall estimation of the method while increasing the computational time by a factor of 4. Also, it is known that the change of variables in Equation (13) present a _gauge freedom_, meaning that any exact gradient of the form \(\nabla\psi\) added to \(\Phi\), will end in the same velocity profile. In this work, we did not explore how

Figure 10: Velocity streamlines of the transient problem for the reference solution (a and d), the mean estimated velocity estimated by the PINNs (b and e), and the velocity obtained after a FEM simulation using the found parameters (c and f) at two time instants: \(t=0.12\ s\), which correspond to peak systole, and at \(t=0.45\ s\), during mid diastole.
different gauges could impact the inverse problem performance. Future work may include studying those and how they can help enhance the properties of the optimization method [56].

The present methodology is not free of limitations and has to be seen as a first step in using neural networks for solving clinical real-world inverse problems. One of the limitations of this work is the lack of experiments with real clinical data. In this line, future work includes the use and assimilation of real patient MRI blood flow images, which we expect will need better strategies in how to deal with such data, as previous work in [22] and in [57] have done to mitigate the presence of artifacts and poor quality of the images.

Finally, since the flexibility in defining a forward and inverse problem within the PINNs workflow, this method could be easily extended to other clinical-relevant tasks such as the denoising and undersampling of MRI blood flow measurements as it was done in [58; 59; 60] with high dependence on a numerical solver and model assumptions. In this line, the PINNs workflow could be seen as a good candidate with less restriction in using contaminated low-resolution data.

## 5 Conclusions

In this work, we present how a physics-informed neural network (PINNs) can be used for estimating hemodynamical information from MRI-like images in a synthetic scenario, where a reduced-order model was coupled with the 3D Navier-Stokes equations to generate the blood flow images. The PINNs were able to solve the coupled system in two regimes: steady and transient flow, using simulated 2D MRI images, a mean pressure curve, the non-slip condition at the wall, and the model equations. Furthermore, at the end of the training, the algorithm was able to estimate the best matching physical parameters, which in this case, represent the elastic response of the adjacent vasculature of the study domain. The overall result shows good accuracy when compared with the reference values in the steady regime and moderately worsened results when dealing with the transient state, where the number of parameters is increased by three-fold, and the data dimensionality is also increased.

To summarize, the proposed framework provides an initial step to truly personalize hemodynamic models from clinical data using a neural network workflow. We envision patient-specific models enabling better diagnostics and treatments for various cardiovascular diseases.

## Acknowledgements

This work was funded by ANID Chile: Millennium Science Initiative Program ICN17_002 (IMFD), ICN2021_004 (iHealth) and NCN19_161 (ACIP), Fondecyt grant 11201250, FONDECYT-Iniciacion 11220816, Basal Funds for Center of Excellence FB210005 (CMM) and Fondecyt Postdoc grant 3230549.