# Parameter estimation in fluid flow models from undersampled frequency space data

Miriam Löcke, Pim van Oojj, Cristóbal Bertoglio

<figure><img src="image_1.png" /></figure>

## To cite this version:

Miriam Löcke, Pim van Oojj, Cristóbal Bertoglio. Parameter estimation in fluid flow models from undersampled frequency space data. 2025. hal-04973326

HAL Id: hal-04973326

https://hal.science/hal-04973326v1

Preprint submitted on 3 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Parameter estimation in fluid flow models from
undersampled frequency space data

Miriam Löcke¹, Pim van Oojj², Cristóbal Bertoglio¹

¹Bernoulli Institute, University of Groningen, The Netherlands

²Department of Radiology and Nuclear Medicine, Amsterdam University Medical Centers,
University of Amsterdam, The Netherlands

March 3, 2025

Abstract

4D Flow MRI is the state of the art technique for measuring blood flow, and it provides valuable information for inverse problems in the cardiovascular system. However, 4D Flow MRI has a very long acquisition time, straining healthcare resources and inconveniencing patients. Due to this, usually only a part of the frequency space is acquired, where then further assumptions need to be made in order to obtain an image.

Inverse problems from 4D Flow MRI data have the potential to compute clinically relevant quantities without the need for invasive procedures, and/or expanding the set of biomarkers for a more accurate diagnosis. However, reconstructing MRI measurements with Compressed Sensing techniques introduces artifacts and inaccuracies, which can compromise the results of the inverse problems. Additionally, there is a high number of different sampling patterns available, and it is often unclear which of them is preferable.

Here, we present a parameter estimation problem directly using highly undersampled frequency space measurements. This problem is numerically solved by a Reduced-Order Unscented Kalman Filter (ROUKF). We show that this results in more accurate parameter estimation for boundary conditions in a synthetic aortic blood flow than using measurements reconstructed with Compressed Sensing.

We also compare different sampling patterns, demonstrating how the quality of the parameter estimation depends on the choice of the sampling pattern. The results show a considerably higher accuracy than an inverse problem using velocity measurements reconstructed via compressed sensing. Finally, we confirm these findings on real MRI data from a mechanical phantom.
# 1 Introduction

In cardiovascular modeling in general and blood flows in particular, the personalization of spatially distributed (i.e. 3D) models is a key step in performing predictive patient-specific simulations [1]. This requires the estimation of relevant parameters from clinical data, which can be used in diagnostics or in the creation of patient-specific predictive simulations such as 3D models of the hemodynamics in the vascular system.

A common technique for acquiring measurements of blood flows is phase-contrast Magnetic resonance Imaging (MRI), or PC-MRI [2], as it is non-invasive and non-ionizing. However, acquiring 4D flow data with MRI requires long acquisition times, which strains clinical resources and inconveniences patients. This applies especially if a high spatial or temporal resolution in the data is required.

MRI applies sequences of magnetic fields and measures the precession of the quantum spin of hydrogen atoms in a strong magnetic field. The resulting measurement equates to the Fourier transformation of the complex magnetization, the phase of which contains information about the velocity. As such, MRI scans in frequency space (also called k-space). One of the features of MRI is that it is possible to select which frequencies of the k-space to measure. Therefore, to reduce the long acquisition times, often only a small part of the k-space is acquired [3].

This partial acquisition introduces artifacts in the reconstructed image that reduce both the image quality. As a result, various Compressed Sensing (CS) techniques have been developed, which aim to reconstruct velocities from highly undersampled data while minimizing artifacts [4]. CS has allowed for more important reductions in the sampling, in the sense that the accuracy of the reconstructed velocity images is higher than with the previously developed sampling and reconstruction strategies [5]. Nonetheless, as we will show later in this article, using CS-reconstructed data for inverse problems compromises the accuracy of the estimated flow and parameters considerably. Secondly, the potential of directly using the k-space measurements in the inverse problem has remain unexplored, while having the potential to avoid the introduction of reconstruction artifacts.

In this paper, we introduce a new technique for estimating blood flow parameters directly from the undersampled k-space MRI data by using an objective function designed for complex MRI measurements. This both avoids artifacts and improves efficiency by skipping the reconstruction/compressed sensing step.

# 2 Theory

## 2.1 MRI measurement and reconstruction in a nutshell

### 2.1.1 Velocity encoding

Let us denote by $u(x, t)$ the component of the velocity $u(x, t)$ in the direction $d$ (fixed in space and in time), and $x \in \mathbb{R}^3$, $t \in \mathbb{R}$ represent the spatial and temporal location within the image space, respectively. The MR images usually analyzed in the clinical context are
actually complex valued, namely the so-called magnetization

$$m(\boldsymbol{x}, t) = M(\boldsymbol{x}, t)e^{i\phi(\boldsymbol{x}, t)} \in \mathbb{C}, \qquad (1)$$

where the magnitude $M > 0$ usually displays the anatomy, and the phase $\phi(\boldsymbol{x}, t) \in (-\pi, \pi]$ in blood flow imaging takes the form

$$\phi(\boldsymbol{x}, t) = \frac{\pi}{venc} u(\boldsymbol{x}, t) + \phi_{back}(\boldsymbol{x}, t) \qquad (2)$$

where $\phi_{back}$ is the background phase, and $venc \in \mathbb{R}$ is the selected velocity encoding (which is a function inversely proportional to the strength and duration of the velocity encoding magnetic gradients).

Since both $u$ and $\phi_{back}$ are unknown on every voxel, at least two measurements are required. In the simplest case (which is the most commonly used in clinical practice), one measurement with no gradient is acquired to obtain $\phi_{back}$, and a second one with $venc \neq 0$, leading to two complex magnetic measurements, so that the velocity can be simply obtained by subtracting the phases, namely,

$$\hat{u} = \frac{\phi - \phi_{back}}{\pi} venc. \qquad (3)$$

Therefore, this technique is called Phase-Contrast MRI (PC-MRI).

In practice, the image $m(\boldsymbol{x}, t)$ is discrete in both space and time, being divided into voxels and time instants. A *voxel* refers to a spatial unit across which the signal is assumed to be constant to provide a single value. The number of voxels is determined by the spatial resolution. The time instants are the times at which a measurement is made, which are determined by the temporal resolution. Therefore $m(\boldsymbol{x}, t)$ can be represented as a matrix in $\mathbb{C}^{N_x \times N_y \times N_z \times N_T}$, with $N_x, N_y, N_z$ being the number of voxels in each spatial direction and $N_T$ the number of time instants. We will refer to the total number of voxels as $N = N_x \cdot N_y \cdot N_z$.

## 2.1.2 The raw signal and image reconstruction

As a matter of fact, the raw signal measured by the MRI scanner corresponds to the spatial Fourier transform of the magnetization, which can, at each measured “time instant”, be formulated as

$$Y(\boldsymbol{k}) = \mathcal{F} \left[ \boldsymbol{M} \odot \exp \left( i \left( \frac{\pi}{venc} \boldsymbol{u}_{meas} + \phi_{back} \right) \right) \right] (\boldsymbol{k}) + \epsilon(\boldsymbol{k}) \qquad (4)$$

where $\odot$ denotes the Hadamard product. We now consider discrete quantities $\boldsymbol{M} \in \mathbb{R}^N$, $\boldsymbol{u}_{meas} \in \mathbb{R}^N$ and $\mathcal{F}: \mathbb{R}^N \to \mathbb{R}^N$ is the three-dimensional discrete Fourier transform defined as

$$\mathcal{F}[\boldsymbol{X}](\boldsymbol{k}) = \sum_{n_x=0}^{N_x-1} \sum_{n_y=0}^{N_y-1} \sum_{n_z=0}^{N_z-1} \boldsymbol{X}_{n_x n_y n_z} \exp \left( -i2\pi \left( \frac{k_1 n_x}{N_x} + \frac{k_2 n_y}{N_y} + \frac{k_3 n_z}{N_z} \right) \right) dx \qquad (5)$$
and $\epsilon(\mathbf{k}) \in \mathbb{C}^N$ is a complex Gaussian noise with a mean of zero. Hence the measurements
are $\mathbf{Y}^n \in \mathbb{C}^N$ for $n = 1, \dots, N_T$.

In practice, for a given temporal resolution, only a number of k-space lines can be ac-
quired. Therefore, the MRI scanner acquires different lines in k-space data at each cardiac
cycles. Therefore a “time instant” *n* actually contains frequencies measured during differ-
ent cardiac cycles, under the assumption that there are no significant differences between
consecutive cardiac cycles (or if there are differences, for instance in the heart rate, such
data is rejected). There is also freedom in choosing different k-space trajectories at different
instants of the cardiac cycle, though that is rarely done.

The noise level depends on the setup chosen in the scans. For instance, the lower the
spatial resolution (i.e., the larger the voxel size), the higher the signal-to-noise ration (because
the signal increases, for the same noise level of the device). Also, the lower the *venc* the
higher the sensitivity of the magnetization phase to the velocity, and therefore one expects
lower noise in the reconstructed velocities with phase-contrast. As the phase is limited to the
range [−π, π), the velocity is however limited to the range [−*venc*, *venc*). If the maximum
velocity exceeds the chosen *venc*, this results in phase aliasing artifacts. As a result, a good
choice of *venc* necessitates a compromise between noise level and presence of artifacts.

If the frequency space is fully sampled, i.e. the values of all k-space locations are known,
the velocity can be reconstructed by applying the inverse Fourier transform to $\mathbf{Y}^n$ and
equation (3). However, fully sampling the k-space is not practicable in cardiovascular MRI
scans since 3D scans may take of the order of hours, depending on the spatiotemporal
resolution chosen.

Therefore, to reduce the length of the scan, the frequency space is generally undersampled.
There are a large number of k-space sampling patterns used in MRI applications, usually in
combination with a reconstruction algorithm specific to the sampling pattern.

Some common 2D patterns include regular Cartesian sampling, radial sampling, and
pseudo-spiral sampling [6]. These can be extended into 3D by stacking 2D designs, leading
to the stack-of-stars and stack-of-spirals designs [6]. The patterns can differ in whether lines
are arranged in a horizontal and vertical or in a shear pattern, in the number of radial
lines or spiral arms, the number of points per line, and whether the patterns are stacked
with equidistant or variable density in the z-direction. Additionally, with stack-of-stars and
stack-of-spirals, there is the option of rotating the pattern by a certain angle for each level
in the stack.

Other methods recommend to use pseudo-random sampling methods, usually with a
higher density in the center, to ensure that the artifacts are incoherent and resemble noise.
Here as well there are many choices for which probability distribution to pick the points
from.

Undersampled data cannot be reconstructed with an inverse Fourier transform, as it
violates the Nyquist limit. *Compressed sensing* (CS) techniques are commonly used to re-
construct undersampled data in a way that limits the presence of artifacts. These techniques
often exploit the sparsity of the data in a different domain, e.g the wavelet domain. Generally
this can be formulated as solving the optimization problem

$$ \hat{\mathbf{m}} = \underset{\mathbf{m}}{\operatorname{argmin}} \{ \| \mathcal{F}_U(\mathbf{m}) - \mathbf{Y} \|_2^2 + \lambda \| \Phi \mathbf{m} \|_1 \} \quad (6) $$

where $\mathcal{F}_U$ is the Fourier transform combined with the undersampling, $\mathbf{Y}$ is the measured k-space data, $\hat{\mathbf{m}}$ is the reconstructed image data, $\Phi$ is a sparsifying transform, and $\lambda$ is the regularization parameter.

The sparsifying transform transports the image into a domain where the desired image is sparse. For MRI images, several such domains exist, for example wavelet domains or total variation time. A good choice of $\Phi$ depends on the spatial and temporal resolution as well as the nature of the artifacts. This makes pseudo-random masks appealing, as they lead to incoherent noise artifacts. The regularization parameter $\lambda$ is usually determined heuristically or empirically and it is known to have an impact in the reconstructed images.

## 2.2 The k-space parameter estimation problem

### 2.2.1 Parameters estimation in PDEs

To consider the underlying physics of the the application, we assume that they can be modelled as the solution $\mathbf{u} : \mathbb{R}^3 \times [0, T] \to \mathbb{R}^3$ of a partial differential equation

$$ F \left( \mathbf{u}, \frac{\partial \mathbf{u}}{\partial x_1}, \frac{\partial \mathbf{u}}{\partial x_2}, \frac{\partial \mathbf{u}}{\partial x_3}, \frac{\partial \mathbf{u}}{\partial t}, \boldsymbol{\theta} \right) = 0 \quad (7) $$

in a domain $\Omega \subset \mathbb{R}^3$ with an initial condition $\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}^0$ and a set of boundary conditions which is dependent on a set of parameters $\boldsymbol{\theta}$, which can describe boundary conditions or material parameters. Then we can define the forward operator $\mathcal{A}(\boldsymbol{\theta})$ which describes the solution of this PDE for a vector of model parameters $\boldsymbol{\theta} \in \mathbb{R}^p$, which includes material parameters, values of boundary conditions, etc. As such, the forward problem generates data according to a physical model with given model parameters.

The goal of the inverse problem overlying this PDE is to retrieve a (sub)set of the parameters $\boldsymbol{\theta}$ given measurements of $\mathbf{u}(\mathbf{x}, t)$.

### 2.2.2 Measurements

For formulating the inverse problem, it can be seen that we have three different choices for the measurements:

1. the frequency space data $\mathbf{Y}^n \in \mathbb{C}^N$ for $n = 1, \dots, N_T$, which is complex-valued, raw data which might be undersampled. The noise on this data is Gaussian and zero-mean.

2. the magnetization data $\mathbf{M}_{\text{meas}}^n \in \mathbb{C}^N$ for $n = 1, \dots, N_T$, which relates to the frequency data as $\mathbf{M}_{\text{meas}}^n = \mathcal{F}^{-1}(\mathbf{Y}^n) = \mathbf{M}_{\text{true}} + \mathcal{F}^{-1}(\boldsymbol{\epsilon})$. This is also complex-valued data with a Gaussian noise with zero mean, as the inverse Fourier transform of a Gaussian is Gaussian as well. However, if the k-space is undersampled, this measurement has to
be reconstructed with CS. In this case, the noise on the velocity measurements cannot assumed to be Gaussian and independent, as there are correlations between the noise in different pixels. [7]

3. the velocity data $V^n \in \mathbb{R}^N$ for $n = 1, \dots, N_T$. $V^n$ is calculated from the angle of the magnetization, meaning that in the case of undersampling the magnetization has to be reconstructed first using compressed sensing as well. Additionally, the noise on this measurement cannot be approximated well by a Gaussian distribution[8].

For all three types of measurements, each velocity direction is acquired separately, therefore it is possible that only one or two velocity components are available rather than all three, or that they are acquired in a direction that does not match a vector component in the canonical basis.

### 2.2.3 Objective function

In a Bayesian framework, the inverse problem for the estimation of parameters from velocity measurements can be solved by minimizing the functional

$$ \hat{\theta} = \underset{\theta \in \mathbb{R}^p}{\operatorname{argmin}} \frac{1}{2\sigma_v^2} \sum_{n=1}^{N_T} \sum_{s=1}^{N} \left( [V^n - \mathcal{H}(\mathbf{u}_{\theta}^n)]_s \right)^2 + \frac{1}{2} \|\theta - \theta^0\|_{(\mathbf{P}^0)^{-1}}^2 \quad (8) $$

where $V^k \in \mathbb{R}^N$ are the measurements of the velocity provided by the reconstruction of the PC-MRI acquisition and $\mathcal{H}: [H_1(\Omega)]^3 \to \mathbb{R}^N$ is the observation operator, which is applied to the result $\mathbf{u}_{\theta}^n$ of the forward model for a set of model parameters $\theta \in \mathbb{R}^p$ at the time corresponding to the measurement $n$. In this standard inverse problem for velocities, this observation operator usually corresponds to the interpolation from the model's (usually finer) mesh to an array of measured velocities at each image voxel. In fact, $[\cdot]_s$ denotes the s-th vector element (or voxel in the case of image data), therefore summing over all the vector elements. $\theta^0$ is the initial guess for the parameters with its covariance matrix $\mathbf{P}^0$, both of which are given by the user but they are assumed to be known for the numerical inverse solver. $\sigma_v$ denotes the standard deviation of the noise on the velocity measurements, which we also assume to be (approximately) known.

This inverse problem has been commonly applied to the estimation of boundary conditions in hemodynamics. For a thorough review, see [1].

In [9], an alternative formulation of this objective function was proposed to account for aliasing artifacts present in velocity MRI data. The functional in this case is

$$ \begin{aligned} \hat{\theta} &= \underset{\theta \in \mathbb{R}^p}{\operatorname{argmin}} \frac{1}{2\sigma_M^2} \sum_{n=1}^{N_T} \sum_{s=1}^{N} \left( [\Re(\mathbf{M}_{\text{meas}}^n) - |\mathbf{M}_{\text{meas}}^n| \cos(\phi_{\text{back}} + \mathcal{H}(\mathbf{u}_{\theta}^n) \frac{\pi}{venc})]_s \right)^2 \\ &\quad + \frac{1}{2\sigma_M^2} \left( [\Im(\mathbf{M}_{\text{meas}}^n) - |\mathbf{M}_{\text{meas}}^n| \sin(\phi_{\text{back}} + \mathcal{H}(\mathbf{u}_{\theta}^n) \frac{\pi}{venc})]_s \right)^2 + \frac{1}{2} \|\theta - \theta^0\|_{(\mathbf{P}^0)^{-1}}^2 \\ &= \underset{\theta \in \mathbb{R}^p}{\operatorname{argmin}} \sum_{n=1}^{N_T} \sum_{s=1}^{N} \frac{|\mathbf{M}_{\text{meas}}^n|_s^2}{\sigma_M^2} \left( 1 - \cos(\frac{\pi}{venc}([V^n - \mathcal{H}(\mathbf{u}_{\theta}^n)]_s)) \right) + \frac{1}{2} \|\theta - \theta^0\|_{(\mathbf{P}^0)^{-1}}^2 \quad (9) \end{aligned} $$
where $M_{meas}^n$ are measurements of the complex magnetization and $\sigma_M$ is the standard deviation of the noise for each of the (complex magnetization) real/imaginary measurements. This formulation assumes that the magnitude $|M_{meas}^n|$ of the magnetization is known in order to formulate the problem in terms of the measured and observed velocities.

Still, the aforementioned cost functions do not account for artifacts originating from the frequency undersampling of the data. Therefore, as we will observe in the numerical examples, the error in the data and hence in the estimated parameters in the inverse problem drastically grows when increasing the undersampling. Here, we postulate to solve this issue by formulating the parameter estimation problem as it is done in CS through a data fidelity term in terms of the k-space, therefore leading to the minimization problem:

$$ \hat{\theta} = \underset{\theta \in \mathbb{R}^p}{\operatorname{argmin}} \frac{1}{2\sigma_y^2} \sum_{n=1}^{N_T} \sum_{s=1}^{N} \left( [\Re(Y^n - \mathcal{H}_F(\mathbf{u}_\theta^n))]_s \right)^2 + \left( [\Im(Y^n - \mathcal{H}_F(\mathbf{u}_\theta^n))]_s \right)^2 + \frac{1}{2} \|\theta - \theta^0\|_{(\mathcal{P}^0)^{-1}}^2 \quad (10) $$

with the observation operator $\mathcal{H}_F$ being defined as

$$ \mathcal{H}_F(u) = \mathcal{F} \left( \mathbf{M} \odot e^{i \frac{\pi}{v_{enc}} \mathbf{u} + \phi_{back}} \right) \odot \mathbf{S} \quad (11) $$

where $\mathbf{S} \in \mathbb{R}^N$ is the sampling mask, i.e. with the entries corresponding to the sampled voxels set to 1, and the others to 0. The $Y^n$ are measurements in frequency space and $\sigma_y$ is the standard deviation of the noise on the frequency-space measurements, which is assumed to be Gaussian and zero-mean. The use of this observation operator requires that the magnitude of the magnetization, or an approximation thereof, is known.

### 2.2.4 Parameter estimation with Kalman filtering

There are different approaches to solving the optimization problems described above. Adjoint-based variational data assimilation approaches fit the entire ensemble of measurements to the ensemble of the observations at the matching time steps. However, this poses significant storage requirements, as the entire trajectory needs to be stored. On the other hand, sequential data assimilation assimilates the information from each new measurement at the point that measurement occurs during one forward pass of the model. This leads to a sequential improvement of the estimation over the duration of the inverse problem as more measurements are considered.

We have chosen to use the Reduced Order Unscented Kalman Filter (ROUKF) in order to solve the optimization problem (10), as it is computationally tractable and has been successfully used in blood-flow problems and other time-resolved problems due to its recursivity.

ROUKF is a sequential parameter estimator which corrects the posterior distribution of the state and parameters in each time step with the available measurement. In order to do so, ROUKF generates a set number of particles sampled from the prior distribution, and propagates each of them through the forward problem. The result of the propagation for each of the particles' state is used to calculate the correction of the mean and the covariance of the parameters at each time measurements are available.
As a result, the filter relies on the measurement error or innovation, which is proportional to the derivative of the data fidelity term of the objective function with respect to $u$. For the objective function (8), this leads to the innovation

$$ \Gamma^n = V^n - \mathcal{H}(u^n) \qquad (12) $$

Similarly, for the objective function in (9), the innovation results from the derivative of the cost function with respect to the state, namely

$$ \Gamma^n = \frac{1}{\sqrt{2}} |\mathcal{M}(t^n)| \sin \left( \frac{\pi}{venc} \cdot (\mathbf{V}^n - \mathcal{H}(\mathbf{u}^n)) \right) $$

where the factor $\frac{\sqrt{2venc}}{\pi} |\mathcal{M}(t^n)|^{-1}$ is added to ensure the equivalence of this innovation with Equation (12) in the case of high $venc$. For more details, see [9].

As done for the first cost function, assuming Gaussianity in complex magnetization measurements, the innovation for the k-space based objective function (10) is defined by:

$$ \Gamma_n = \begin{bmatrix} \Re(\mathbf{Y}^n) - \Re(\mathcal{H}_F(\mathbf{u}^n)) \\ \Im(\mathbf{Y}^n) - \Im(\mathcal{H}_F(\mathbf{u}^n)) \end{bmatrix}. \qquad (13) $$

Moreover, we remark that we are adapting the Kalman filter as described in [10]. Instead of simplex sigma points (which involve $p+1$ solutions of the forward problem), we are using canonical sigma points, which involve $2p$ solutions of the forward problem but they have shown in our simulations better performance than their simplex counterparts. In particular, the results using both types of sigma points deviate more the larger the number of parameters to be estimated. Also, the results with simplex points change depending on the enumeration of the parameters during the estimation process, which is not the case with canonical points. Additionally, to ensure positivity of the parameters, we are reparameterizing such that $\theta = \theta^0 2^{\nu}$ with $\theta^0$ the initial guess for the parameters and the filter being applied to $\nu$, as is also done in [9].

# 3 Methods

## 3.1 Synthetic data

### 3.1.1 Forward problem setup

We use the same model as is used in [9]. We consider a geometry of the lumen of the ascending and descending aorta including the outlets of the brachiocephalic artery, left common carotid artery, and left subclavian artery, as depicted in Figure 1. This geometry serves as the domain for the forward model.

The boundary of the geometry consists of six different boundaries: $\Gamma_{in}$ being the inlet boundary in the ascending aorta, $\Gamma_w$ the arterial wall, and the remaining boundaries $\Gamma_l$ for $l = 1, \dots, 4$ representing the outlets.
<figure><img src="image_3.png" /><figcaption>Figure 1: 3D aortic model geometry</figcaption></figure>

We model the blood flow in this domain with the incompressible Navier-Stokes equations for the velocity $\mathbf{u}(\mathbf{x}, t) \in \mathbb{R}^3$ and the pressure $p(\mathbf{x}, t) \in \mathbb{R}$:

$$ 
\left\{
\begin{array}{l}
\rho \frac{\partial \mathbf{u}}{\partial t} + \rho(\mathbf{u} \cdot \nabla)\mathbf{u} - \mu \Delta \mathbf{u} + \nabla p = 0 \text{ in } \Omega \\
\nabla \cdot \mathbf{u} = 0 \text{ in } \Omega \\
\mathbf{u} = \mathbf{u}_{in} \text{ on } \Gamma_{in} \\
\mathbf{u} = \mathbf{0} \text{ on } \Gamma_{w} \\
\mu \frac{\partial \mathbf{u}}{\partial n} - p \mathbf{n} = -P_l(t) \mathbf{n} \text{ on } \Gamma_l
\end{array}
\right.
\quad (14) $$

with $\rho, \mu$ the density and dynamic viscosity of the fluid and $P_l(t)$ being given by a Windkessel boundary condition defined by:

$$ 
\left\{
\begin{array}{l}
P_l = R_{p,l}Q_l + \pi_l \\
Q_l = \int_{\Gamma_l} \mathbf{u} \cdot \mathbf{n} dx \\
C_{d,l} \frac{d\pi_l}{dt} + \frac{\pi_l}{R_{d,l}} = Q_l
\end{array}
\right.
\quad (15) $$

This boundary condition models the effects of the remaining vascular system on the outlet via the proximal and distal resistances $R_p, R_d$ of the vasculature and the distal compliance $C_d$ of the vessels.

The inflow $\mathbf{u}_{in}$ is defined as

$$ \mathbf{u}_{in} = -U f(t) \mathbf{n} $$

where $U$ is a constant amplitude and

$$ f(t) = \begin{cases} \sin(\frac{\pi t}{T}) & \text{if } t \le T \\ \frac{\pi}{T}(t-T) \exp^{-k(t-T)} & \text{if } T_c > t > T \end{cases} $$
with $T_c = 0.8$ and $T = 0.36$.

The physical parameters are set as seen in Table 1. The forward problem is solved using

<table><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>ρ (gr·cm<sup>3</sup>)</td><td>1.2</td></tr><tr><td>μ (P)</td><td>0.035</td></tr><tr><td>U (cm·s<sup>-1</sup>)</td><td>75</td></tr><tr><td>T<sub>c</sub> (s)</td><td>0.80</td></tr><tr><td>T (s)</td><td>0.36</td></tr><tr><td>κ (s<sup>-1</sup>)</td><td>70</td></tr></tbody></table>

<table><caption>Table 1: Physical parameters and numerical values of the three-element Windkessel parameters for every outlet.</caption><thead><tr><th></th><th>Γ<sub>1</sub></th><th>Γ<sub>2</sub></th><th>Γ<sub>3</sub></th><th>Γ<sub>4</sub></th></tr></thead><tbody><tr><td>R<sub>p</sub> (dyn·s·cm<sup>-5</sup>)</td><td>480</td><td>520</td><td>520</td><td>200</td></tr><tr><td>R<sub>d</sub> (dyn·s·cm<sup>-5</sup>)</td><td>7200</td><td>11520</td><td>11520</td><td>4800</td></tr><tr><td>C (dyn<sup>-1</sup>·cm<sup>5</sup>)</td><td>4·10<sup>-4</sup></td><td>3·10<sup>-4</sup></td><td>3·10<sup>-4</sup></td><td>4·10<sup>-4</sup></td></tr></tbody></table>

a semi-implicit 3D-0D coupling scheme as seen in [9]. The full algorithm is detailed in the appendix.

### 3.1.2 Synthetic measurements

The forward solution is generated with a time step of $dt = 1$ms and undersampled in time to $dt_{meas} = 15$ms, leading to a total of 56 measurements. From the solution of the forward problem, we simulate a PC-MRI acquisition by subsampling into a rectangular measurement mesh with a resolution of [2mm, 2mm, 2mm] and then applying the process described in Section 2.1 with a *venc* of double the maximal velocity. The magnitude is modelled as

$$M(\boldsymbol{x}) = \begin{cases} 1.0 & \text{if } \boldsymbol{x} \text{ is in the lumen of the vessel} \\ 0.5 & \text{otherwise.} \end{cases} \qquad (16)$$

Finally a complex Gaussian noise $\epsilon \in \mathbb{C}^N$ is added with a signal-to-noise ratio (SNR) of 15. Fifty independent realizations of the noise were generated.

For comparison, we reconstructed velocity measurements from these synthetic measurements using the Berkeley Advanced Reconstruction Toolbox (BART)[11]. BART is a command-line-based software that provides a flexible framework of compressed sensing methods, as well as tools for simulation, pre-processing, and image reconstruction, providing a multitude of different regularization options. In this work, we have used this toolbox for compressed sensing reconstructions of the velocities, using total variation in time as for the regularization.

Next, the sampling mask is applied to these simulated frequency space measurements. We take a 2D subsampled mask in the $x-y$-plane and sample fully in the $z$-direction as in [4]. We consider different subsampling rates $R = \frac{N_{sampled}}{N_{total}} = 8, 16, 32$, with two different masks: the pseudo-spiral mask and the pseudo-random Gaussian mask, which is sampled according to a Gaussian probability distribution, as shown in Figure 2. For the pseudo-spiral mask, the points are placed evenly on a cartesian grid along a spiral with six turns and a final radius reaching the edge of the mask.
Additionally, we require measurements of the magnitude itself.

As anatomical images are usually readily available and the magnitude generally does not depend on the encoding direction, making a cheap 2D/3D acquisition feasible, we consider a reconstructed magnitude from the measurements with $R = 2$ and a gaussian mask using temporal l1-regularization in BART with a regularization parameter $\lambda = 0.001$. Examples of the phase and magnitude of the measurements are shown in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 2: Sampling masks</figcaption></figure>

### 3.1.3 Inverse problem setup

Without pressure measurements, not all the Windkessel resistances in the system can be uniquely determined at once. Therefore we have fixed the values for $\Gamma_4$ (the outlet in the descending aorta), and estimate the amplitude of the inflow $U$ as well as the distal resistances $R_{d,k}$, $k = 1, 2, 3$ of the remaining Windkessel boundary conditions.

We consider two different initial guesses for these parameters, as in [9]

1. "low" guess: $U = 40, R_{d,1} = 4000, R_{d,2} = 4000, R_{d,3} = 4000$

2. "high" guess: $U = 150, R_{d,1} = 20000, R_{d,2} = 20000, R_{d,3} = 20000$

while the target values are $U = 75, R_{d,1} = 7200, R_{d,2} = 11520, R_{d,3} = 11520$.

The initial standard deviation for the reparameterized parameters was set to 0.5, i.e. $P^0 =$
<figure><img src="image_5.png" /><figcaption>(a) K-space phase, fully sampled (b) True magnitude of the magnetization</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>(c) True velocity in z-direction</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) K-space magnitude, fully sampled (e) Magnitude of the magnetization reconstructed with BART (f) Velocity in z-direction reconstructed with BART from R=2 from R=1</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 3: Examples of simulated measurements, taken at a slice in the z-direction.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Gaussian mask, $R = 8$, $\lambda = 0.01$ (b) Gaussian mask, $R = 16$, $\lambda = 0.01$ (c) Gaussian mask, $R = 32$, $\lambda = 1.0$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Spiral mask, $R = 8$, $\lambda = 0.1$ (e) Spiral mask, $R = 16$, $\lambda = 1.0$ (f) Spiral mask, $R = 32$, $\lambda = 10.0$</figcaption></figure>

Figure 4: Examples of velocities reconstructed with BART using different masks and acceleration factors. Depending on the sampling mask, different kinds of artifacts appear in the reconstructed velocity.
0.5I, meaning that the prior models that there is ≈ 95% probability that the target value will lie between the range half/twice the initial guess. The standard deviation of the noise was estimated from the initial time step by computing the standard deviation of

$$Y^0 - M(t^0) \odot \exp (i\phi_{back}(t^0)) \odot S$$

which is the noise of the data under the assumption that the velocity at the first time step is zero. The estimated standard deviation for each case is shown in Table 2. The estimation of the noise in this case is impacted by noise in the reconstructed magnitude.

<table><caption>Table 2: Estimated standard deviation of the noise for each mask and acceleration factor, using the magnitude reconstructed from R2 for acceleration factors greater than 1</caption><thead><tr><th>Acceleration factor</th><th>Gaussian</th><th>Spiral</th></tr></thead><tbody><tr><td>R1/actual value</td><td colspan="2">15.526</td></tr><tr><td>R8</td><td>15.829</td><td>13.813</td></tr><tr><td>R16</td><td>18.716</td><td>13.791</td></tr><tr><td>R32</td><td>23.043</td><td>13.844</td></tr></tbody></table>

## 3.2 Phantom data

### 3.2.1 Flow phantom measurements

We used data from a flow experiment on a phantom of the carotid artery reported in [4]. The phantom is made of a distensible silicon and is suspended in water. It takes the shape of a bifurcating tube, to simulate a carotid artery, and a backflow tube. A pump generates a pulsatile flow with a rate of roughy 60 bpm, simulating a cardiac cycle lasting 1s.

The 4D Flow MRI scan was performed with a 3T Ingenia scanner by Philips Healthcare. All three velocity directions were acquired, plus one acquisition with no encoding gradient to acquire the background phase. The scan parameters were set to TR = 8.9ms, TE = 4.5ms, FA = 8°, venc=150cm/s. The matrix size was [160, 160, 40] with a spatial resolution of [0.8mm, 0.8mm, 0.8mm] and a temporal resolution of 19 frames per cardiac cycle. The scan was accelerated with an acceleration factor of R=2. The scan used a 32-channel coil leading to 15 independent coil measurements, resulting in a total matrix size of [160, 160, 40, 4, 15, 19].

### 3.2.2 Forward problem setup

A structured mesh was created from a segmentation of the lumen, matching the spatial resolution of the scan. From this an expanded, higher-resolution mesh was created using Blender, Meshlab and gmsh. This mesh was used for the forward problem. The two meshes are depicted in Figure 6.

We again model the blood flow with the incompressible Navier-Stokes equations, with the same choices for the physical parameters. The inflow boundary is modelled by a Dirichlet
<figure><img src="image_13.png" /><figcaption>(a) K-space magnitude of (b) x-component of the (c) Zero-filled magnitude of (d) Zero-filled magnitude of coil 4 velocity reconstructed with the magnetization of coil 4 the magnetization of coil 6 BART from R=2</figcaption></figure>

Figure 5: Examples of the phantom measurements. Reconstructed magnitudes of the magnetization of different coils show the different coil sensitivities.

boundary condition with an inflow

$$ \boldsymbol{u}_{in} = \boldsymbol{u}_{profile} U f(t) $$

where $\boldsymbol{u}_{profile}$ is a flow profile provided by the solution of a Stokes problem in the domain, $U$ is a constant amplitude, and $f(t)$ defines the temporal profile as

$$ f(t) = \begin{cases} \sin(\frac{\pi t}{T}) & \text{if } t \le \frac{3}{4}T \\ \sin(\frac{3}{4}\pi)(1 - t + \frac{3}{4}T)\exp\left(-\left(t - \frac{3}{4}T\beta\right)\right) & \text{if } t > \frac{3}{4}T \end{cases} $$

with $T = 0.64$, $\beta = 5$. The shape was determined empirically to match the shape of the flow rate in the phantom.

The two outlets were modelled with resistance boundary conditions

$$ P_l = R_{p,l}Q_l \text{ on } \Gamma_l $$

$$ Q_l = \int_{\Gamma_l} \boldsymbol{u} \cdot \boldsymbol{n} dx $$

As the flow division between the outlets is determined by the proportion of the resistances, we can fix $R_{p,1} = 100$ at an arbitrary value.

### 3.2.3 Inverse problem setup

We estimate the inflow amplitude $U$ and the resistance $R_{p,2}$ of the boundary condition at the right outlet. The initial guess is $U = 100$, $R_{p,2} = 100$, i.e. an equal division of flow between the outlets. The initial standard deviation was set to 0.5.
<figure><img src="image_14.png" /><figcaption>Figure 6: Meshes of the carotid phantom</figcaption></figure>
We are reparameterizing such that $\theta = \theta^0 \odot \nu$ with $\theta^0$ the initial guess for the parameters and the filter being applied to $\nu$. This was done because the exponential reparameterization used above proved unstable in this case.

To extend the filter to multiple coils, each of the coil measurements was treated as a separate, independent measurement, resulting in 15 measurements per time step. The variance of the noise for each coil was estimated from the first measurement by computing the standard deviation of

$$Y^0 - M(t^0) \odot \exp(i\phi_{back}(t^0)) \odot S$$

as before, which neglects potential spatial variation of the variance due to the sensitivity of the coils.

The observation operator requires measurements of the magnitude of the magnetization and the background phase. The background phase is reconstructed as

$$\phi_{back} = \angle \mathcal{F}^{-1}(Y_{back}) \qquad (17)$$

where $Y_{back}$ is a zero-filled measurement acquired with $R = 2$ and no encoding gradient, i.e. capturing only the background magnetization and no fluid velocity.

For the magnitude, we consider two options. The first option is

$$M = |\mathcal{F}^{-1}(Y)| \qquad (18)$$

for each velocity direction, where $Y$ is a zero-filled measurement acquired with $R = 2$ in that velocity direction. This however assumes that a highly sampled measurement is available. Therefore the other option is to use

$$M = |\mathcal{F}^{-1}(Y_{back})| \qquad (19)$$

for all velocity directions, using the same measurement as for $\phi_{back}$. This would only require a highly sampled measurement of one out of four encoding gradients.

The original data with an acceleration factor of 2 were acquired with an incoherent pseudo-spiral sampling in the $k_y - k_z$ direction with a fully sampled $k_x$-direction. We further undersample this by applying masks with a Gaussian or spiral sampling pattern onto the sampling mask for $R = 2$ to achieve higher acceleration factors. In contrast to the masks for the synthetic data, this results in different masks for each time step and velocity direction. Examples of the resulting masks are shown in Figure 7.

The measurements are assumed to be placed at every 0.053s. For the inverse problem, the first two measurements (at 0.053s and 0.106s) are omitted as the presence of negative velocity values in these measurements could interfere with the Kalman filter, since this could lead to negative values in the particles which would be unphysical for the forward simulation.
<figure><img src="image_15.png" /><figcaption>(a) Original, $R = 2$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Spiral, $R = 8$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Spiral, $R = 16$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Spiral, $R = 32$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(e) Spiral, $R = 64$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(f) Spiral, $R = 128$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(g) Gaussian, $R = 8$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(h) Gaussian, $R = 16$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(i) Gaussian, $R = 32$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(j) Gaussian, $R = 64$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(k) Gaussian, $R = 128$</figcaption></figure>

Figure 7: Sampling masks for the phantom data, taken at the $k_x = 0$ slice at time step 0
# 4 Results

## 4.1 Synthetic data

In order to accurately reflect the different impact of the parameters on the flow, we compute the error in terms of the flow reconstructed from the estimated parameters. Using the estimated parameters, we solve the forward problem 14 again with a time step of $dt = 0.001$. The error is calculated as

$$e = \frac{\|\mathbf{u}_{ref} - \mathbf{u}_{recon}\|_2}{\|\mathbf{u}_{ref}\|_2}$$

where $\mathbf{u}_{recon}$ and $\mathbf{u}_{ref}$ are vectors consisting of the velocity values for all components at each point in the geometry at each point in time ($dt_{meas} = 15ms$) stacked together, using the estimated parameter values for $\mathbf{u}_{recon}$ and the true parameter values for $\mathbf{u}_{ref}$.

We compare the results of our method (directly estimating from frequency-space measurements) to velocity measurements reconstructed with compressed sensing using BART with a 1-regularization in time. The regularization parameter was determined empirically by visual evaluation for each combination of acceleration factor and mask and is listed in Table 3. For the estimation, we are using the method presented in [9] to take advantage of knowing the $R = 2$ magnitude.

<table><caption>Table 3: Regularization values $\lambda$ used for the reconstruction with BART.</caption><thead><tr><th>Acc. Factor</th><th>Spiral</th><th>Gaussian</th></tr></thead><tbody><tr><td>R1</td><td></td><td>0.0001</td></tr><tr><td>R8</td><td>0.1</td><td>0.01</td></tr><tr><td>R16</td><td>1.0</td><td>0.01</td></tr><tr><td>R32</td><td>10.0</td><td>1.0</td></tr></tbody></table>

The error values for all different acceleration factors for the low and high initial guess are depicted in Figure 8, for our method as well as using BART, as well as the error values of the reconstructed velocity itself (without first solving the inverse problems). In all cases, using the inverse problem to retrieve parameters and then reconstruct the flow achieves better results than the velocity measurements themselves. For the k-space cost function, the errors of both masks are very close to the error with fully sampled data for $R = 8$, and increase for $R = 16$ and $R = 32$. The gaussian mask achieves lower error values for all three acceleration factors.

Considering the velocity measurements from data reconstructed using BART, it can be seen that the gaussian mask performs considerably better than the spiral mask, especially for $R = 8$. This matches the expectation, as a pseudo-random masks leads to incoherent artifacts, which can be easily excluded with a temporal regularizer. Nonetheless, the error increases drastically for $R = 16$ and $R = 32$. The spiral mask shows a strong increase in error for $R = 8$ already and remains high with higher subsampling rates.

Using frequency measurements directly outperforms using the reconstructed velocity measurements for all acceleration factors except for the fully sampled case.
<figure><img src="image_26.png" /><figcaption>Figure 8: Error values for different acceleration factors and different masks. Dotted lines are the error values of the CS-reconstructed velocity values themselves, dashed lines for the flow reconstructed from the inverse problem from BART velocities, solid lines are from the inverse problem in frequency space. Low initial guess on the left, high initial guess on the right. The bars indicate the standard deviation of the error.</figcaption></figure>

Both the high and the low initial guess show the same pattern. As such, from here on, we will limit ourselves to showing results for the low initial guess as the ones for the high initial guess provide no additional information.

By comparing the estimated parameter values of each parameter individually, as seen in Figure 9, it is apparent that the success of the estimation differs depending on the parameter. The inflow $U$ is estimated relatively accurately by all masks, with only a small decrease with increasing $R$. In comparison, the distal resistance of the first Windkessel outlet ($R_{d,1}$) is underestimated considerably and increasingly by the spiral mask, while the gaussian mask estimates it well for $R = 8$ and $R = 16$, but overestimates it severely for $R = 32$. $R_{d,2}$ and $R_{d,3}$ show the same pattern, with the spiral mask underestimating the values compared to the gaussian mask, though more severely for $R_{d,2}$ than $R_{d,3}$. For $R = 32$, the spiral mask slightly overestimates the value of $R_{d,3}$.

We also consider the case where measurements are only available for one of the velocity components, in this case the $z$-component, which is equivalent to the foot-head direction in this setup. Figure 10 shows the error values for this case. The error values show the same patterns as when using all velocity components, and are only slightly higher. Here there is a lesser difference in the error values between the spiral and gaussian mask for $R = 8$ and $R = 16$, and the spiral mask shows a flatter error curve, whereas the error for the gaussian mask increases with $R = 32$. Again, using the frequency measurements achieves lower error values than BART measurements for all subsampled data.

When considering each parameter separately in Figure 11, the same patterns and similar
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>Figure 9: Estimated values of individual parameters. The dashed line indicates the true value.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 10: Error values for different acceleration factors using only the *z*-component of the velocity.</figcaption></figure>

values persist as for all velocity components. The exception is the Gaussian mask for $R_{d,2}$, which now significantly underestimates the values, unlike the spiral mask.

## Robustness to the choice of *venc*

As described previously in Remark 2, a low *venc* leads to a high signal-to-noise ratio, but can lead to aliasing artifacts if the actual maximal velocity exceeds the *venc*. The cost function used in [9] remedies this by distinguishing the actual velocity from the wrapped ones when the physical parameters affect the velocity on several voxels simultaneously. In the present work, the cost function corresponds to a similar cost function as in [9] – but with an additional Fourier transform – therefore also including the aliasing compensation.

To investigate the robustness of our method to choosing *venc*s lower than the maximal velocity, thus utilizing the higher sensitivity to the velocity in the signal, we compare three additional *venc* values corresponding to 80%, 30%, and 10% of the maximum velocity. The results are shown in Figure 12.

The results improve considerably with a decreasing *venc* for 80% and 30%, especially for higher acceleration factors. However, the smallest *venc* of 10% leads to increased error values and standard deviations, with the inverse problem even failing to provide results for the spiral mask for $R = 32$. The results for both masks perform similarly with the decreasing *venc*, indicating that this response is inherent to the inverse problem formulation – as shown
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>Figure 11: Estimated values of individual parameters using only the *z*-component of the velocity. The dashed line indicates the true value.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>(a) Spiral</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(b) Gaussian</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(c) BART with spiral mask</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(d) BART with gaussian mask</figcaption></figure>

Figure 12: Error values for different values of the *venc*. Spiral mask on the left, Gaussian mask on the right. The value for *R32* for the spiral mask for the lowest *venc* is excluded as the inverse problem failed to converge. For the BART reconstructions, *venc* = 10% was not considered.
in [9] – rather than a result of the choice of mask.

In comparison, with the BART reconstructions, the error remains similar for *venc* = 80%*V*<sub>max</sub> and then increases significantly for *venc* = 30%*V*<sub>max</sub>. This behaviour can be explained from the fact that the reconstructed velocities were already unwrapped prior to the parameter estimation as described in [9].

Using an estimated magnitude of the magnetization

In the prior results, we have been using a magnitude reconstructed from k-space subsampled
with *R* = 2. To check the accuracy of this substitution, we now compare to providing perfect
knowledge of the magnitude. We also compare to using a constant magnitude with a value
of 0.5, to model having little-to-no knowledge of the magnitude.

<figure><img src="image_40.png" /><figcaption>Figure 13: Error values for different estimations of the magnitude</figcaption></figure>

As can be seen in Figure 13, using the magnitude reconstructed from frequency data with an acceleration factor of *R* = 2 achieves results that are very close to using the perfect magnitude. On the other hand, using a constant magnitude leads to a significant increase in error that appears constant across the different acceleration factors. In this case, there are also less differences between the two masks. Nonetheless, for higher acceleration factors, even with a constant magnitude the error is less than with velocity measurements reconstructed from BART.

We also provide examples of the parameter curves of ROUKF in Figure 14, which show
the development of the estimated parameter over the length of the simulation. The shaded
area indicates the variance of the parameter. This parameter variance, for each of the masks
using our new method, is higher at the end of the simulation for higher *R*, since less data
is available. This effect is not observed for the BART-reconstructed velocities, where the
variance decreases similarly for all acceleration factors. Here the variance is also much smaller
<figure><img src="image_41.png" /><figcaption>(a) Spiral, $U$</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(b) Gaussian, $U$</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(c) BART (Gaussian mask), $U$</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>(d) Spiral, $R_{d,1}$</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(e) Gaussian, $R_{d,1}$</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>(f) BART (Gaussian mask), $R_{d,1}$</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>(g) Spiral, $R_{d,2}$</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>(h) Gaussian, $R_{d,2}$</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(i) BART (Gaussian mask), $R_{d,2}$</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(j) Spiral, $R_{d,3}$</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(k) Gaussian, $R_{d,3}$</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(l) BART (Gaussian mask), $R_{d,3}$</figcaption></figure>

Figure 14: Evolution of the estimated parameters over the time of the simulation. The dashed line indicates the true value.
than for the frequency-space based method, indicating higher confidence in the estimated
parameter. This higher confidence is based on the fact that due to the reconstruction, the
measurements contain more data than in frequency space (as they now have as many voxels
as the image space, as compared to only the actually measured voxels). Since these data
however originate from the same subsampled measurements, they do not actually contain
more information, which leads to this misleading quantification of the uncertainty in the
inverse problem.

Furthermore, it can be seen that for low $R$, the curves qualitatively follow the shape
of the curve for the fully sampled data, but with higher sampling rates, the choice of the
mask not only quantitatively but also qualitatively impacts the evolution of the parameter
as it is estimated. The curves from our method also differ considerably from those using
BART, which develop erratic spikes in either direction as the acceleration factor increases.
In comparison to that, the curves for each mask remain more similar to each other with
increasing $R$.

It can also be seen that in all the cases, the value of the estimated parameter is mostly
constant after roughly $t = 0.5$, and the parameter variance does not decrease considerably
after that. This indicates that for both frequency-space data and reconstructed data, the
relevant data is located entirely within systole.

## 4.2 Phantom data

The results of the estimation of the noise per mask and acceleration factor are shown in Figure 15. While the estimation of the noise with the spiral mask remains close to constant for the acceleration factors, the estimations with the gaussian mask show major variations for some of the coils. This is likely due to the spatial variation of the signal strength for some of the coils and the fact that the gaussian mask does not uniformly cover the space. Nonetheless this approach provides a good approximation of the standard deviation of the noise even for high acceleration factors.

We are using the same error metric as described in Section 4.1, with $u_{ref}$ the solution of a forward problem with the parameters estimated from $R = 2$.

The error values can be seen in Figure 16. A considerable difference between the behaviour of the two masks is apparent. Both achieve similar results for an acceleration factor of $R = 16$, but the error increases higher for the spiral mask than for the Gaussian mask. The Gaussian mask also shows a decrease in error from $R = 64$ to $R = 128$, which may be a result of differences in the noise estimation, or otherwise an outlier. Both choices for which magnitude to use in the observation show similar results, implying robustness to using a magnitude from a different velocity direction.

In Figure 17, we show the flow rates of the reconstructed flows at the inlet and each of the outlets. It can be seen that both masks tend to overestimate the inflow with increasing acceleration factor $R$. For the gaussian mask, the flow split between the two outlets remains fairly consistent as $R$ increases, whereas the spiral mask no longer accurately depicts the flow split from $R = 32$ upwards. This can also be observed by looking at the estimated parameter values in Figure 18. Both masks show a very consistent estimation of the inflow
<figure><img src="image_53.png" /><figcaption>(a) Gaussian mask</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(b) Spiral mask</figcaption></figure>

Figure 15: Estimation of the standard deviation of the signal noise per coil, for various acceleration factors

<figure><img src="image_55.png" /></figure>

Figure 16: Error values for different acceleration factors for both error metrics. The "mag0" labelled results use the magnitude from the acquisition without a velocity encoding gradient.
amplitude *U*, with differences only for *R* = 128, but the estimation of the resistance boundary condition *R*<sub>*p*,2</sub> shows larger changes for different acceleration factors. The spiral mask especially overestimates *R*<sub>*p*,2</sub> for higher acceleration factors, leading to the lack of distinction of the outlets in the flow rates.

<figure><img src="image_56.png" /><figcaption>Figure 17: Flow rates at the inlet and each of the outlets of the forward simulations based on the estimated parameters.</figcaption></figure>

Examples of the ROUKF parameter curves are provided in Figure 19. It can be seen
that the curves for the gaussian and spiral masks have a similar shape with some qualitative
differences. The acceleration factor seems to have a higher effect on the parameter evolution
<figure><img src="image_57.png" /><figcaption>Figure 18: Estimated parameter values for different *R* for each mask. Errorbars depict the parameter variance, multiplied by a factor of 1000 for visibility. The "mag0" labelled results use the magnitude from the acquisition without a velocity encoding gradient.</figcaption></figure>

curve for the spiral mask than the gaussian mask. It can also be seen that the parameter variance decreases drastically with the first measurement, with only a very small variance afterward despite continuing changes in the parameter values.

# 5 Discussion

The numerical results are consistent between the simulated dataset and the phantom dataset. In both cases, our approach robustly estimates boundary conditions of the flow even with high undersampling rates and low *venc* values and outperforms the conventional approach of reconstructing velocity measurements from undersampled data through Compressed Sensing.

The results also show significant differences between the two studied sampling patterns, with qualitatively and quantitatively different outcomes of the inverse problem. The sampling patterns also appear to differ in their estimation of different parameters, indicating that the choice of sampling pattern may depend on which parameter should be estimated more accurately.

A limitation of the method is that it requires an accurate estimate of the magnitude and the background phase of the magnetization. However, since the magnitude can be scanned without motion encoding gradients, a highly sampled image of the magnitude is quick to acquire. We have also shown that it was possible – at least in the phantom dataset – to acquire the background image highly sampled for approximations of the magnitude and background phase, while undersampling the remaining three directions. As we have shown some robustness of the method with regards to flaws in the magnitude, one could also acquire the magnitude only at diastole, and use this for all time steps to cut down on acquisition
<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /><figcaption>Figure 19: Evolution of the estimated parameters over the time of the simulation. Crosses mark times where corrections according to the measurements are made.</figcaption></figure>
time, or correct it using the model itself. Additionally, the noise present in the reconstructed
magnitude and background phase is currently not considered by the inverse problem. This
will be subject of future research.

Another limitation of this paper is that real patient data is not included, which would
pose additional complexities such as the movement of the vessels due to the breathing of the
patient.

# 6 Conclusion

We have proposed a new formulation for the inverse problem for parameter estimation in
fluid flow problems using undersampled frequency-space MRI data and demonstrated it
using a Reduced-Order Unscented Kalman Filter. This method outperforms the results of
an inverse problem using velocity data reconstructed with Compressed Sensing, especially
for high undersampling rates and different *venc* values.

The choice of subsampling mask is shown to have a strong influence on the result of the
estimation of some parameters. Future work could therefore address the challenge of finding
optimal sampling masks for certain parameters.

# Acknowledgments

We acknowledge the funding from the European Research Council (ERC) under the Euro-
pean Union's Horizon 2020 research and innovation program (grant agreement No 852544 -
CardioZoom).

# A Numerical solution method of the forward problem

Here we detail the algorithm used to solve the incompressible Navier-Stokes equation with
Windkessel boundary conditions for the forward problem.
**Algorithm 1:** Fractional step algorithm with a modified semi-implicit Windkessel model coupling

Given the initial conditions $\mathbf{u}^0 = \mathbf{u}(0) \in V_{\Gamma_{w,h}}$ and $\pi_1^0, \dots, \pi_N^0 \in \mathbb{R}$, perform for $j > 0$, with $t^j = j\tau$:

1. **Viscous Step:** Find the tentative velocity $\tilde{\mathbf{u}}^n \in V_{\Gamma_{w,h}}$ such that:

$$ \left\{ \begin{aligned} & \tilde{\mathbf{u}}^j|_{\Gamma_{in}} = \mathbf{u}_{inlet}(t^j) \\ & \frac{\tau}{\rho}(\tilde{\mathbf{u}}^j, \mathbf{v})_{\Omega_h} + \rho(\mathbf{u}^{j-1} \cdot \nabla \tilde{\mathbf{u}}^j, \mathbf{v})_{\Omega_h} + \frac{\rho}{2}((\nabla \cdot \mathbf{u}^{j-1})\tilde{\mathbf{u}}^j, \mathbf{v})_{\Omega_h} + (\delta\mathbf{u}^{j-1} \cdot \nabla \tilde{\mathbf{u}}^j, \mathbf{u}^{j-1} \cdot \nabla \mathbf{v})_{\Omega_h} \\ & +2\mu(\epsilon(\tilde{\mathbf{u}}^j), \epsilon(\mathbf{v}))_{\Omega_h} + \sum_{\ell=1}^{K} \frac{\rho}{2}|\mathbf{u}^{j-1} \cdot \mathbf{n}|_{-}(\tilde{\mathbf{u}}^j, \mathbf{v})_{\Gamma_{\ell}} = \frac{\tau}{\rho}(\mathbf{u}^{j-1}, \mathbf{v})_{\Omega_h} \end{aligned} \right. \qquad (20) $$

for all $\mathbf{v} \in V_{\Gamma_{in} \cup \Gamma_{w,h}}$, and $|x|_{-}$ denoting the negative part of $x$.

2. **Projection-Windkessel Step:** Compute $\tilde{Q}^j = \int_{\Gamma_{\ell}} \tilde{\mathbf{u}}^j \cdot \mathbf{j}$. Find $p^j \in Q_h$ such that:

$$ \frac{\tau}{\rho}(\nabla p^j, \nabla q)_{\Omega_h} + \sum_{\ell=1}^{K} \frac{\overline{p^j}_{\Gamma_{\ell}} \overline{q}_{\Gamma_{\ell}}}{\gamma_{\ell}} + \epsilon \sum_{\ell=1}^{K} (\mathcal{T}(\nabla p^j), \mathcal{T}(\nabla q))_{\Gamma_{\ell}} = \sum_{\ell=1}^{K} \left( \tilde{Q}^j + \frac{\alpha_{\ell} \pi_{\ell}^{j-1}}{\gamma_{\ell}} \right) \overline{q}_{\Gamma_{\ell}} - (\nabla \cdot \tilde{\mathbf{u}}^j, q)_{\Omega_h}, \qquad (21) $$

for all $q \in Q_h$ and with $\overline{(\cdot)}_{\Gamma_{\ell}} = \frac{1}{\text{Area}(\Gamma_{\ell})} \int_{\Gamma_{\ell}} (\cdot) ds$ and $\mathcal{T}(\mathbf{f}) = \mathbf{f} - (\mathbf{f} \cdot \mathbf{j})\mathbf{j}$.

3. **Velocity correction Step:** Find $\mathbf{u}^j \in [L^2(\Omega_h)]^3$ such that:

$$ (\mathbf{u}^j, \mathbf{v})_{\Omega_h} = (\tilde{\mathbf{u}}^j - \frac{\tau}{\rho} \nabla p^j, \mathbf{v})_{\Omega_h} $$

for all $\mathbf{v} \in [L^2(\Omega_h)]^3$

4. **Update-Windkessel Step:** Set $P_{\ell}^j = \overline{p^j}_{\gamma_{\ell}}$ and compute $\pi_{\ell}^j \in \mathbb{R}$ as:

$$ \pi_{\ell}^{j} = \left( \alpha_{\ell} - \frac{\alpha_{\ell} \beta_{\ell}}{\gamma_{\ell}} \right) \pi_{\ell}^{j-1} + \frac{\beta_{\ell}}{\gamma_{\ell}} P_{\ell}^{j}, \quad \ell = 1, \dots, K $$

## B ROUKF algorithm

Here we detail the ROUKF algorithm adapted from [10]. Let us first consider the notation $[Z^{(*)}]$ as the matrix with the column-wise collection of vectors $Z_{(1)}, Z_{(2)}, \dots$

Define the *canonical sigma-points* $I_{(1)}, \dots, I_{(2p)} \in \mathbb{R}^p$ such that

$$ I_{(i)} = \begin{cases} \sqrt{p}e_i, & \text{for } 1 \le i \le p \\ -\sqrt{p}e_{i-p}, & \text{for } p+1 \le i \le 2p \end{cases} $$

where the vectors $e_i$ form the canonical base of $\mathbb{R}^p$. Moreover, define the weight $\alpha = \frac{1}{2p}$.
We denote by $\hat{\mathbf{X}}_-, \hat{\mathbf{X}}_+ \in \mathbb{R}^r$ a priori (model prediction) and a posteriori (corrected by observations) estimates of the true state $\mathbf{X}^n \in \mathbb{R}^r$. In the semi-implicit coupled 3D-0D fractional step Algorithm 1, the state consists in the velocity field $\mathbf{u}^n$ and the Windkessel pressures $\pi^n$. Estimates of all unknown parameters are summarized by the corresponding a priori and a posteriori vectors $\hat{\theta}_-, \hat{\theta}_+ \in \mathbb{R}^p$. The discretized forward model is written as $\mathbf{X}^n = \mathcal{A}^n(\mathbf{X}^{n-1}, \theta^{n-1}), \mathcal{A}^n$ denoting the model operator.

For given values of the initial condition $\hat{\mathbf{X}}_+^0 = \mathbf{X}^0 \in \mathbb{R}^r$, the initial expected value of the parameters $\hat{\theta}_+^0 = \theta^0 \in \mathbb{R}^p$ and its covariance matrix $\mathbf{P}^0$, perform

* Initialization: initialize the sensitivities as

$$ \mathbf{L}_{\theta}^{0} = \sqrt{\mathbf{P}^{0}} \quad (\text{Cholesky factor}), \quad \mathbf{L}_{X}^{0} = \mathbf{0} \in \mathbb{R}^{r \times p}, \quad \mathbf{U}_{0} = \mathbf{P}_{\alpha} \equiv \alpha [\mathbf{I}^{(*)}][\mathbf{I}^{(*)}]^{T} \quad (22a) $$

Then, for $n = 1, \dots, N_T$:

* Sampling: generate $2p$ particles from the current state and parameter estimates, i.e. for $i = 1, \dots, 2p$:

$$ \left\{ \begin{aligned} \hat{\mathbf{X}}_{(i)}^{n-1} &= \hat{\mathbf{X}}_{+}^{n-1} + \mathbf{L}_{X}^{n-1}(\mathbf{C}^{n-1})^{T}\mathbf{I}_{(i)}, \\ \hat{\theta}_{(i)}^{n-1} &= \hat{\theta}_{+}^{n-1} + \mathbf{L}_{\theta}^{n-1}(\mathbf{C}^{n-1})^{T}\mathbf{I}_{(i)} \end{aligned} \right. \quad (22b) $$

with $\mathbf{C}^{n-1}$ the Cholesky factor of $(\mathbf{U}^{n-1})^{-1}$.

* Prediction: propagate each particle with the forward model and compute an a priori state prediction:

$$ \left\{ \begin{aligned} \hat{\mathbf{X}}_{(i)}^{n} &= \mathcal{A}^{n}(\hat{\mathbf{X}}_{(i)}^{n-1}, \hat{\theta}_{(i)}^{n-1}), \quad \hat{\theta}_{(i)}^{n} = \hat{\theta}_{(i)}^{n-1}, \quad i = 1, \dots, 2p \\ \hat{\mathbf{X}}_{-}^{n} &= E_{\alpha}([\hat{\mathbf{X}}^{n}^{(*)}]) \equiv \alpha \sum_{i=1}^{2p} \hat{\mathbf{X}}_{(i)}^{n} \\ \hat{\theta}_{-}^{n} &= E_{\alpha}([\hat{\theta}^{n}^{(*)}]) \end{aligned} \right. \quad (22c) $$

* Correction: compute a posteriori estimates based on measurements for state and parameters, using the $i$-th particle innovation $\Gamma_{(i)}^n$:

$$ \left\{ \begin{aligned} \mathbf{L}_{X}^{n} &= \alpha[(\hat{\mathbf{X}}^{n}^{(*)})[\mathbf{I}^{(*)}]^{T} \\ \mathbf{L}_{\theta}^{n} &= \alpha[(\hat{\theta}^{n}^{(*)})[\mathbf{I}^{(*)}]^{T} \\ \mathbf{L}_{\Gamma}^{n} &= \alpha[(\Gamma^{n}^{(*)})[\mathbf{I}^{(*)}]^{T} \\ \mathbf{U}^{n} &= \mathbf{P}_{\alpha} + (\mathbf{L}_{\Gamma}^{n})^{T}\mathbf{W}^{-1}\mathbf{L}_{\Gamma}^{n}, \quad \mathbf{P}_{\alpha} = \alpha[\mathbf{I}^{(*)}][\mathbf{I}^{(*)}]^{T} \\ \hat{\mathbf{X}}_{+}^{n} &= \hat{\mathbf{X}}_{-}^{n} - \mathbf{L}_{X}^{n}(\mathbf{U}^{n})^{-1}(\mathbf{L}_{\Gamma}^{n})^{T}\mathbf{W}^{-1}E_{\alpha}([\Gamma^{n}^{(*)}]) \\ \hat{\theta}_{+}^{n} &= \hat{\theta}_{-}^{n} - \mathbf{L}_{\theta}^{n}(\mathbf{U}^{n})^{-1}(\mathbf{L}_{\Gamma}^{n})^{T}\mathbf{W}^{-1}E_{\alpha}([\Gamma^{n}^{(*)}]) \\ \mathbf{P}_{\theta}^{n} &= \mathbf{L}_{\theta}^{n}(\mathbf{U}^{n})^{-1}(\mathbf{L}_{\theta}^{n})^{T} \\ \mathbf{P}_{X}^{n} &= \mathbf{L}_{X}^{n}(\mathbf{U}^{n})^{-1}(\mathbf{L}_{X}^{n})^{T} \end{aligned} \right. \quad (22d) $$

with $\mathbf{W} = \sigma_y^2 \mathbf{I}$.
References

[1] David Nolte and Cristóbal Bertoglio. Inverse problems in blood flow modeling: A review. *International journal for numerical methods in biomedical engineering*, 38(8):e3613, 2022.

[2] Michael Markl, Frandics P. Chan, Marcus T. Alley, Kris L. Wedding, Mary T. Draney, Chris J. Elkins, David W. Parker, Ryan Wicker, Charles A. Taylor, Robert J. Herfkens, and Norbert J. Pelc. Time-resolved three-dimensional phase-contrast mri. *J Magn Reson Imaging*, 17(4):499–506, Apr 2003.

[3] Lukas M Gottwald, Eva S Peper, Qinwei Zhang, Bram F Coolen, Gustav J Strijkers, Aart J Nederveen, and Pim van Ooij. Pseudo-spiral sampling and compressed sensing reconstruction provides flexibility of temporal resolution in accelerated aortic 4d flow mri: A comparison with k-t principal component analysis. *NMR in Biomedicine*, 33(4):e4255, 2020.

[4] Eva S Peper, Lukas M Gottwald, Qinwei Zhang, Bram F Coolen, Pim van Ooij, Aart J Nederveen, and Gustav J Strijkers. Highly accelerated 4d flow cardiovascular magnetic resonance using a pseudo-spiral cartesian acquisition and compressed sensing reconstruction for carotid flow and wall shear stress. *Journal of Cardiovascular Magnetic Resonance*, 22(1):1–15, 2020.

[5] Michael Lustig, David Donoho, and John M. Pauly. Sparse mri: The application of compressed sensing for rapid mr imaging. *Magnetic Resonance in Medicine*, 58(6):1182–1195, 2007.

[6] Li Feng. Golden-angle radial mri: Basics, advances, and applications. *Journal of Magnetic Resonance Imaging*, 56(1):45–62, 2022.

[7] Lauren Partin, Daniele E. Schiavazzi, and Carlos A. Sing Long. An analysis of reconstruction noise from undersampled 4d flow mri, 2022.

[8] Pablo Irarrazaval, Ali Dehghan Firoozabadi, Sergio Uribe, Cristian Tejos, and Carlos Sing-Long. Noise estimation for the velocity in mri phase-contrast. *Magnetic Resonance Imaging*, 63:250–257, 2019.

[9] Jeremías Garay, David Nolte, Miriam Löcke, and Cristóbal Bertoglio. Parameter estimation in fluid flow models from aliased velocity measurements. *Inverse Problems*, 38(9):095002, aug 2022.

[10] P. Moireau and D. Chapelle. Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems. *COCV*, 17:380–405, 2011. doi:10.1051/cocv/2010006.
[11] Moritz Blumenthal, Martin Heide, Christian Holme, Martin Juschitz, Bernhard Rapp, Philip Schaten, Nick Scholand, Jon Tamir, Christian Tönnes, and Martin Uecker. mrire-con/bart: version 0.9.00, December 2023.