Journal Pre-proof Transient flow prediction in an idealized aneurysm geometry using data assimilation Franziska Gaidzik, Daniel Stucht, Christoph Roloff, Oliver Speck, Dominique Thévenin, Gábor Janiga
PII: DOI: Reference:
S0010-4825(19)30371-3 https://doi.org/10.1016/j.compbiomed.2019.103507 CBM 103507
To appear in:
Computers in Biology and Medicine
Received date : 9 July 2019 Revised date : 27 September 2019 Accepted date : 12 October 2019 Please cite this article as: F. Gaidzik, D. Stucht, C. Roloff et al., Transient flow prediction in an idealized aneurysm geometry using data assimilation, Computers in Biology and Medicine (2019), doi: https://doi.org/10.1016/j.compbiomed.2019.103507. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Conflict of Interest Statement
Journal Pre-proof Conflict of interest
urn al P
Highlights (for review)
Computers in Biology and Medicine
“Transient flow prediction in an idealized aneurysm geometry using data assimilation“
Improved flow prediction in intracranial aneurysms by merging PC-MRI data with CFD simulations
Transient data assimilation using a localization approach with Ensemble Kalman Filtering to reduce ensemble size requirements
Two different update processes for the data assimilation experiment are evaluated
urn al P
Unique validation of PC-MRI and DA with high-resolution transient PIV data to generate ground-truth flow data
4D Flow; Data assimilation; Ensemble Kalman Filter; intracranial aneurysm; PC-MRI
Journal Pre-proof *Revised manuscript (clean) Click here to download Revised manuscript (clean): Revised_Clean.tex
Transient flow prediction in an idealized aneurysm geometry using data assimilation
Franziska Gaidzika , Daniel Stuchtb,c , Christoph Roloffa , Oliver Speckb,d , Dominique Th´evenina , G´abor Janigaa,∗ a
Lab. of Fluid Dynamics and Technical Flows, University of Magdeburg “Otto von Guericke” b Institute of Experimental Physics, University of Magdeburg “Otto von Guericke” c Institute of Biometry and Medical Informatics, University of Magdeburg “Otto von Guericke” d Leibniz Institute for Neurobiology, Magdeburg, Germany
Hemodynamic simulations are restricted by modeling assumptions and uncertain initial and boundary conditions, whereas Phase-Contrast Magnetic Resonance Imaging (PC-MRI) data is affected by measurement noise and artifacts. To overcome the limitations of both tech-
urn al P
niques, the current study uses a Localization Ensemble Transform Kalman Filter (LETKF) to fully incorporate noisy, low-resolution Phase-Contrast MRI data into an ensemble of highresolution numerical simulations. The analysis output provides an improved state estimate of the three-dimensional blood flow field in an intracranial aneurysm model. Benchmark measurements are carried out in a silicone phantom model of an idealized aneurysm under pulsatile inflow conditions. Validation is ensured with high-resolution Particle Imaging Velocimetry (PIV) obtained in the symmetry plane of the same geometry. Two data assimilation approaches are introduced, which differ in their way to propagate the ensemble members in time. In both cases the velocity noise is significantly reduced over the whole cardiac cycle. Quantitative and qualitative results indicate an improvement of the flow field prediction in
comparison to the raw measurement data. Although biased measurement data reveal a systematic deviation from the truth, the LETKF is able to account for stochastically distributed errors. Through the implementation of the data assimilation step, physical constraints are ∗
Corresponding author Email address: [email protected]
, Universit¨ atsplatz 2, Magdeburg, D-39106, Germany (G´abor Janiga)
Preprint submitted to Computers in Biology and Medicine
September 27, 2019
introduced into the raw measurement data. The resulting, realistic high-resolution flow field can be readily used to asses further patient-specific parameters in addition to the velocity distribution, such as wall shear stress or pressure.
Keywords: 4D Flow, Data assimilation, Ensemble Kalman Filter, intracranial aneurysm,
CFD – Computational Fluid Dynamics; DA – Data assimilation; EnKF – Ensemble Kalman Filters; ICP – Iterative Closest Point; LETKF – Localization Ensemble Transform Kalman Filter; PC-MRI – Phase-Contrast Magnetic Resonance Imaging; PIV – Particle
Image Velocimetry; SNR – signal-to-noise-ratio
Rupture of intracranial aneurysms can cause severe outcomes, such as irreversible dis-
abilities or even death. To improve our understanding of cardiovascular diseases, hemody-
namic studies are very helpful. Such a gain of knowledge can support physicians concerning
outcome prediction and therapy planning. Flow-dependent parameters, such as wall shear
stress or the oscillatory shear index, have been proposed as relevant quantities related to
the growth and rupture probability of aneurysms [1, 2, 3]. Different modalities are available
to assess patient-specific hemodynamic flow fields. Numerical investigations using Compu-
tational Fluid Dynamics (CFD) have been frequently used to identify key parameters and
flow characteristics. However, such simulations require accurate initial and boundary condi-
tions, and depend on the assumptions used in the numerical model [4, 5, 6]. Although high
temporal and spatial resolution can be reached, uncertain initial and boundary conditions
lead to a limited clinical acceptance of the simulation results [7, 8, 9]. Phase-Contrast Mag-
netic Resonance Imaging (PC-MRI) in-vivo measures blood flow quantitatively by encoding
the velocities in the phase of the acquired MR signal [10, 11, 12]. However, in addition to
measurement noise and artifacts, the currently achievable temporal and spatial resolution
urn al P
limits the extraction of clinically-relevant flow features from PC-MRI measurement data.
Alternative measurement techniques, such as Particle Image Velocimetry (PIV), can acquire
high resolution flow fields in-vitro, but a clinical application would be limited or even im-
possible . Being an optically-based technique, PIV requires fully transparent settings.
Patient-specific silicone models allow PIV measurements but are expensive in terms of manu-
facturing time and costs. Additionally, the surface of patient-based silicone models is usually
extracted from MR images, which limits the accuracy of the geometry bounding the flow.
This, together with the pulsatile nature of patient vessels, leads to differences between the
physiological condition and the flow measured within silicone models. Therefore, highly re-
solved PIV can be mainly used as a validation procedure for specific cases and benchmark
Uncertainties associated with either numerical or experimental approaches naturally sug-
gest the use of data assimilation to incorporate the inherently noisy measurement data of
limited resolution into high-resolution numerical simulations. In this manner, accuracy and
physical correctness of the measured velocity field can be improved. De-noising techniques,
as well as divergence-free filtering approaches, can also improve the physical correctness
of the velocity field; however, spatial and temporal resolution remain unchanged [14, 15].
Hence, data assimilation seems to be a promising remedy to improve resolution while enforc-
ing constraints, such as incompressibility and conservation laws. Contrary to most existing
approaches that use available PC-MRI measurement data only as initial conditions for CFD
simulations [4, 16], data assimilation uses a complete measured 4D flow field to be incor-
porated into numerical simulations. In this manner, the reliability of the intracranial state
estimate together, with the prediction of clinically relevant parameters such as wall shear
stress (WSS), should be improved.
urn al P
For fluid dynamics applications regarding compressible flows promising results have been
achieved with adjoint-based data assimilation [17, 18, 19], even if the resulting computational
costs are in the order of 300 CPUh for one iteration loop. That is because solving the
adjoint equations needs approximately twice as much time as solving the direct equation.
More than 30 loops are often required to reach the desired convergence criteria . D’Elia
et al. [20, 21, 22] and Funke et al.  applied variational data assimilation approaches in
intracranial aneurysms. As an alternative to the sequential technique they minimized the
error between observations of a reference flow and a numerical estimation in terms of a
cost function. The authors stated that the need for linear and adjoint models increases the
computation time by a factor of 50 to 100 in comparison to one single CFD simulation. As
a consequence of these high computational costs, most of the stated numerical studies on
variational data assimilation in intracranial aneurysms currently only address steady-state
flows and/or 2D geometries [17, 18, 19, 20, 21, 22]. Funke et al.  investigated transient
3D flow fields, but with a reduced spatial resolution to maintain acceptable computational
costs. Further limitations of previous hemodynamic DA studies result from 1) synthetically
generated measurement data (replaced by real PC-MRI data in the present study), and/or 2)
the lack of a ground truth, in particular for real 4D flow fields (delivered by PIV experiments
in this work).
urn al P
Although promising results have been achieved with this technique for other fluid dynam-
ical applications [24, 25], little attention has been paid to the sequential Kalman Filters in
intracranial aneurysm modeling. Bakhshinejad et al.  implemented an Extended Kalman
Filter for a pulsatile cardiac flow. However, the lack of localization requires a large amount of
ensemble members to ensure filter convergence. The resulting, high computing times reveal
the need for a sequential data assimilation technique that can converge with a limited amount
of ensemble simulations. A possible approach relies on Proper Orthogonal Decomposition
(POD) to reduce the needed ensemble members ; POD can also be used to quantitatively
analyze PC-MRI flow data . As an alternative, the present study predicts an improved
flow field for measured transient intracranial aneurysm flow using a Localization Ensemble
Kalman Filter (LETKF). In general, Ensemble Kalman Filters (EnKF) sample the system
state and covariance matrices by an ensemble of model states. Thus, the covariance ma-
trices are not calculated directly, but estimated through an ensemble and replaced by the
sample covariance. The uncertainty of the numerical model (background) is estimated and
the Ensemble Kalman Filter can be seen as a Monte-Carlo approximation of the original
Kalman Filter . Usually, a large ensemble size is required to sample the system state in
a statistically representative manner. With the implementation of a localization procedure
in the Ensemble Kalman Filter, the amount of needed ensemble members can be reduced
significantly. Building on top of previous investigations , the current study relies on
transient high-resolution PIV data to provide a ground truth, allowing comparisons with
PC-MRI and CFD data and – even more important – a validation of the developed data
2. Material and Methods
urn al P
(a) numerical model
(c) cardiac cycle trajectory
Figure 1: (a) Surface model of the idealized aneurysm used for phantom manufacturing and CFD discretization; (b) MRI setup with the gear pump (1), wave guide through RF-shield (2), flow meter (3), MR-scanner (4) and the 32-channel head coil with the phantom model (5); (c) pulsatile flow profile provided by the programmable pump .
PC-MRI and PIV data are obtained from a two-component silicone model (Wacker RT
601, Burghausen, Germany) of an idealized intracranial aneurysm, with proximal and distal
vessels of 4 mm. To enable a fully developed flow profile, a long (length 40 cm, diame-
ter 4mm), straight, rigid pipe was installed upstream of the aneurysm inlet. Exactly the
same geometry is used for the calculation of the forward model, i.e., the ensemble of CFD
simulations. To mimic realistic flow conditions, a blood substitute is used in all experiments.
It matches approximately the fluid dynamic properties of real blood (ρ = 1222 kg/m3 ;
η = 4.03 mPa s), and simultaneously the refractive index (=1.4122 at 22◦ C) of the silicone
block, which ensures optimal optical access for PIV measurements. The pulsatile flow, rep-
resenting a cardiac cycle (Figure 1 (c)), is provided by a programmable micro-gear pump
(HNP Mikrosysteme, Schwerin, Germany). To prevent any difference concerning the flow,
the experimental set-up was held completely identical (same construction, pipe length, fit-
tings...) for both PIV and PC-MRI measurements. This allows a proper comparison and
cross-validation between the different modalities. Additionally, the resulting data was regis-
tered with the implementation of an Iterative Closest Point (ICP) algorithm, which enables
rigid body transformation. Difficulties in the registration process arose due to geometric
distortions in the PC-MRI data. These artifacts increase with growing distance from the
measurement center, which was chosen in the center of the aneurysm sack. The distorted
parts of distal and proximal vessels are removed prior to the registration step, to ensure
proper matching of the velocity fields in the center of the aneurysm, as acquired from the
Figure 2 shows the flowchart that is used to get an improved state estimate of the intracra-
nial velocity field. A detailed description of the main steps regarding the overall workflow is
given in the following sections.
2.1. Phase-Contrast Magnetic Resonance Imaging (PC-MRI)
urn al P
The 4D flow data was acquired on a 7 Tesla whole-body MRI system (Siemens Health-
ineers, Forchheim, Germany) in a 32-channel head coil (Nova Medical, Wilmington, MA)
using 4D phase-contrast magnetic resonance imaging (PC-MRI). Hereby, the acquisition se-
quence is based on an RF-spoiled gradient echo with quantitative flow encoding in all three
spatial dimensions [30, 31]. The data was acquired with an imaging matrix of 224x224 px,
a 128x128 mm FoV and a slice thickness of 0.57 mm resulting in an isotropic voxel size of
0.57 mm. The total scan time for 11 temporal phases was approx. 109 minutes. The data
of one time step is acquired in the interval between two subsequent time points and results
in a temporal resolution of 82.8 ms. To correct for eddy currents a reference measurement
was acquired, without activated pump and thus without flow inside the aneurysm phantom
model. The phase of the reference data was subtracted from the phase of the flow data to
urn al P
Figure 2: Technical flow chart of the current data assimilation technique. The proposed workflow is conducted for every time point at which a phase-contrast MR image is available. For time-independent data assimilation the workflow can be performed for all 11 time steps individually; whereas for time-dependent data assimilation the ensemble simulations of the subsequent time step are initialized by the analysis outcome of the previous time point.
obtain purely flow-related phase differences. As the flow information is encoded in the phase
of the complex MR-Signal, a velocity encode parameter (venc) specifies the highest velocity,
uniquely encoded in one phase. For the acquired data, this parameter was set to 1.5 m/s.
The signal-to-noise-ratio (SNR) of the acquired images was calculated using the mean of the
signal inside the aneurysm and the signal density of the background noise and was found to
be SNR ≈ 90. The data was post processed using MeVisLab 2.3.1 (MeVis, Bremen, Ger-
many) and the automated tool described in . This includes noise masking, antialiasing
and conversion to the format of the commercial software package EnSight (ANSYS Inc.,
Canonsburg, PA, USA). The acquired flow data serve as an input for the data assimilation
step in the form of observation data.
2.2. Particle Image Velocitmetry (PIV)
The generation of a quantitative gold standard is one of the main challenges in the
formulation of a suitable data assimilation experiment. Particle Image Velocitmetry (PIV) is
currently the most established flow field measurement technique allowing for highly resolved
flow field data. In this study, a high-speed stereoscopic PIV system was used to obtain
time-resolved flow data at the sagittal plane of the idealized aneurysm geometry. By that, a
ground truth is formulated, which is used to validate the data assimilation experiment. The
entire fluidic setup used for the PC-MRI measurements was used without any modification
as the basis for the PIV study. A thin laser light sheet illuminated small resin tracer particles
doped with Rhodamine B (diameter d = 10.46 ± 0.18 µm, density ρ = 1510 kg/m3 ), which
were given to the fluid. Two high-speed cameras recorded double-frame images (inter-frame
delay of 200 µs) of the particles’ fluorescent light at 500 Hz recording frequency. Hence, after
processing (Davis software 8.4.1, LaVision, Goettingen, Germany), a velocity field snapshot
was available every 20 ms during the entire cardiac cycle. A total of 36 independent cardiac
cycles were covered by the PIV recordings; the final phase-averaged velocity field was used
as validation data for this study.
2.3. Background model
urn al P
In the current study, CFD simulations are used as the numerical background model; it
is needed as an input for the formulation of the data assimilation algorithm. To calculate
the background model, an ensemble of CFD simulations is propagated forward in time to
sample the system state and to provide an estimate of the error covariance matrix related
with the numerical model. A block-structured hexahedral mesh is first generated from the
underlying geometry using ANSYS IcemCFD (ANSYS Inc., Canonsburg, PA, USA) result-
ing in approx. 750,000 finite-volume cells. The results of a mesh-independence analysis did
not show any noticeable changes in flow parameters, so that numerical stability with respect
to the mesh size is ensured. Ensemble simulations are carried out using the open source soft-
ware OpenFOAM 4.0 (OpenCFD Ltd., Bracknell, UK). Blood is treated as an isothermal,
incompressible fluid, with density and viscosity prescribed as in the experimental measure-
ments. The vessel walls are assumed to be rigid and no-slip boundary conditions as well
as a zero-pressure outlet are implemented. To ensure computational stability, the time-
step of the background needs to be smaller than the time step provided from the PC-MRI
measurements; it has been chosen to be 1 ms as suggested by .
2.4. Data Assimilation
urn al P
(a) Time-independent DA
(b) Time-dependent DA
Figure 3: For an easier visualization, only 10 ensemble members are shown: (a) Ensemble trajectories for the time-independent data assimilation. Here, the analysis states are independently calculated at the measurement times, indicated by the red points; (b) Schematic illustration of the time-dependent data
By using data assimilation, numerical model data and PC-MRI measurement data are
merged together for an improved state estimate of the flow field in the intracranial aneurysm.
The inflow condition throughout the cardiac cycle, as well as the three-dimensional ve-
locity field, serve as the uncertain quantities in the present data assimilation experiment.
This provides an interesting alternative to approaches frequently used in past studies, in 9
which PC-MRI data were only used to extract inlet and outlet flow rates applied as initial
and/or boundary conditions in a single CFD calculation. The Local Ensemble Transform
Kalman Filter (LETKF) applied in this paper was originally introduced by Harlim and Hunt
[34, 35] in the field of meteorology. It combines the localization method of the Local En-
semble Kalman Filter (LEKF) of Ott et al.  and the Ensemble Transform Kalman Filter
(ETKF) of Bishop et al. . The analysis ensemble is formed as a weighted average of
the background ensemble mean and the observations. Using background and observation
uncertainties, the weights are determined in such a way that the analysis ensemble mean
fits in an optimal manner the given background and observation probability distributions.
By implementing localization, purely local analyses at each model grid point are obtained,
so that only observations within a local region surrounding the grid point are taken into ac-
count. The localization scheme enables efficient parallel computations of the analysis model
state; additionally, it limits the number of needed ensemble members. Beyond model (CFD
simulations) and noisy measurement (PC-MRI) data, following parameters are used as input
for the DA study:
urn al P
• A localization radius of 3 mm, taking only observations into account when they are found within this radius.
• An observational operator H to map the state variables from the m-dimensional simu-
lation space to the s-dimensional observation space. The current observation operator
H is a spatial binning operator, which downsamples the ensemble velocity vectors to
the MRI grid resolution.
• An s × s dimensional observation error covariance matrix R based on the noise char√ 2 venc acteristics of the measurement data: σv = . π SN R
With an randomly large ensemble size, the error covariance matrix of the underlying numer-
ical model can be sampled correctly. Nevertheless, this results in unacceptably high com-
putational costs when running a large amount of multiple CFD resolutions. By introducing
the localization procedure of the LETKF, a smaller ensemble size can be used to sample 10
the subdomains and to determine local weights, which are used to calculate the analysis.
Using 25 ensemble members the filter was able to converge throughout the whole cardiac
cycle. To calculate an improved state estimate, two different data assimilation experiments
are performed, denoted in what follows time-independent (or Case 1), resp. time-dependent
– with cardiac modeling (or Case 2):
1. Case 1: From the raw PC-MRI data the inflow rate at the 11 measurement time steps
can be extracted. Based on the noise characteristics of the PC-MRI measurements
the reconstructed flow rates at every time point are perturbed to generate 25 differ-
ent inflows. Cubic interpolation randomly connects one of the 25 flow rates of the 11
different time points and calculates a trajectory in between the time steps to gener-
ate the ensemble members. Thus, the resulting CFD ensemble simulations differ in
the underlying inflow trajectory throughout the cardiac cycle. Each trajectory is in-
dividually propagated forward in time throughout the cardiac cycle (Figure 3a). For
every measurement time point all ensemble member states are stored separately and
11 different time-independent analyses are calculated. In this case, the analysis is a
function of the ensemble inflow conditions and the measurement data, but is indepen-
dent from the analysis state of a previous time step. No forecast step in the classical
sense can be performed, because no background model for the cardiac cycle trajectory
urn al P
2. Case 2: Similar to Case 1, the ensemble members for the analysis of the first time
point are generated using the noise characteristics of the measurement data. After the
calculation of the first analysis, the cardiac cycle trajectory in between the current and
next time step is reconstructed using a forecast model in the form of Xt = f (Xt−1 ) +
t . Currently, the simple forward model is a linear function, which serves as initial
condition for the ensemble CFD simulations. Thus, each analysis ensemble velocity
field is propagated forward in time, using the previously calculated analysis velocity
field and a linear inflow trajectory as initial conditions.
2.5. Validation The high-resolution PIV measurements of the center slice inside the aneurysm are set as
a ground truth criteria and used to validate the calculated analyses. Quantitative validation
and comparison is based on the interpolation of velocity values obtained by PIV, PC-MRI, or
data assimilation, respectively, onto a random point cloud distributed over the center slice.
The Root-Mean Square Error (RMSE) between the different modalities is calculated and
normalized (NRMSE) to enable a better comparison of the values throughout the cardiac
v u m u1 X t RMSE = (vAnalysis − vPIV )2 m i=1 NRMSE =
RMSE ¯ PIV v
The Similarity Index (SI) measures the similarity of two independent vector sets and
takes deviations of the angle between both vectors (ASI), as well as the differences in the
normalized magnitudes (MSI) into account. For validation purposes the SI, together with
ASI and MSI, between PIV and either PC-MRI or analyses flow fields is evaluated. Equation
urn al P
(2) shows for example the SI calculation between PIV and analysis: |vAnalysis | vAnalysis · vPIV |v | PIV SI = ASI · MSI = 0.5 · 1 + · 1 − − |vAnalysis | · |vPIV | max (|vAnalysis |) max (|vPIV |) (2) 3. Results
3.1. Qualitative flow field validation
First, for a qualitative flow visualization the velocity distribution along the center slice
inside the aneurysm is shown at every time point for which measurement data are available
(Figure 4). PC-MRI and calculated analyses flow fields are compared with the PIV mea-
surements for either time-independent or time-dependent data assimilation. Overall, a good
agreement can be observed between the analyses flow fields and the measured ground truth. 12
urn al P
Figure 4: Qualitative comparison of the transient flow field prediction for both data assimilation experiments together with the PC-MRI and PIV data. The velocity magnitude of a cross-section cut through the center of the aneurysm model is displayed.
However, the velocity magnitude seems to be underestimated by PC-MRI and analysis in
comparison to the PIV data, especially in the high velocity areas. Additionally, the analysis
velocity field calculated in Case 1 seems to be smoother than the flow field in Case 2, in which
a simplified cardiac model was used to propagate the ensemble members forward in time.
In both cases, the LETKF de-noises the PC-MRI data and increases the spatial resolution
of the three-dimensional flow field. The resulting spatial resolution matches the number of
grid points from the underlying ensemble CFD simulations. Thus, the walls of the analyses
cross-section cuts appear to be smoother and better resolved compared to the raw measure-
ment data of PC-MRI. Thanks to this improvement in spatial resolution and extraction of
a smooth surface boundary, clinically relevant parameters like the Wall Shear Stress (WSS)
can be computed. Figure 5 shows the calculated WSS at five different time-points based
on the time-independent, three-dimensional analysis flow field. An area of high shear stress
at the transition from the aneurysm sack to the distal vessel can be clearly identified at all
Figure 5: WSS distribution based on the time-independent three-dimensional analysis flow field. The spatial
3.2. Quantitative flow field validation
urn al P
resolution matches the resolution from the ensemble CFD simulations.
Figure 6: Velocity magnitudes of PIV are plotted against PC-MRI and analysis magnitudes, respectively.
The data is displayed in a 2D histogram density plot.
The flow fields obtained from the different modalities are interpolated onto a randomly
generated point cloud to enable a quantitative comparison. Measurement time point 4
(velocity peak) has been chosen to exemplarily show the performance of the data assimilation
steps, as well as the accuracy of the PC-MRI measurement. The velocity magnitudes of 14
urn al P
Figure 7: Histograms displaying the relative frequency of the similarity index SI (top row), magnitude similarity index MSI (middle row) and angular similarity index ASI (bottom row) over all compared data points for the different modalities. The first quartile (Q1 ) is marked with a vertical blue line.
the different modalities are first plotted against each other (Figure 6). For an illustrative
quantification, a two-dimensional histogram grid is provided, which is additionally color-
coded by the amount of data points corresponding to a specific location; the brighter the
grid cell, the more data points lie there. The red line indicates a perfect match between both
compared modalities. Regarding the agreement between PIV and PC-MRI, a significant
Figure 8: NRMSE (a) and Q1 (b) values of the three different modalities throughout the cardiac cycle.
difference concerning the scatter around the red line is observed in the low and high velocity
ranges. For velocity magnitudes below 0.25 m/s, many points are found and appear to
be equally distributed around the red line (σ = 0.014); the scatter represents the noise of
the measurements. At higher velocities the points become more and more shifted toward
the area above the red line; thus, the MRI measurement increasingly underestimates the
velocities. About 98% of the velocity magnitudes higher than 0.5 m/s are underestimated
by PC-MRI. A similar result is observed for the velocities calculated in the data assimilation
experiments. Although the spread around the red line and thus the noise of the data is
reduced in comparison to PC-MRI (σ = 0.01), the deviation from the ground truth again
increases at higher velocities. For the DA analyses only 83% (Case 1) and 86% (Case 2) of the
velocity values higher than 0.5 m/s are underestimated. Figure 7 shows the relative frequency
of the similarity index (SI), angular similarity index (ASI) and magnitude similarity (MSI)
found at the random sample points. All histogram plots contain a line representing the Q1
value, which quantifies the point that separates the lower 25% of the frequency distribution
from the higher 75%. The indices are calculated between analysis or PC-MRI data compared
to the PIV measurements. For the PC-MRI data, the SI is higher than 0.75 in 75% of the
data points. A much better agreement is found for DA, with a Q1 equal to 0.9 (Case 1) or
0.88 (Case 2), respectively. A more detailed insight in the comparison of the vector fields
can be achieved by splitting the SI into its angular (Figure 7, middle row) and magnitude
urn al P
part (Figure 7, bottom row), respectively. For both indices the Q1 value is increased thanks
to DA. The MSI is higher than 0.96 in 75% of the data points for both analyses, while the
corresponding value is only 0.92 for the PC-MRI data. This effect is even more pronounced
for the angular differences. While the Q1 value for the PC-MRI data is relatively low (0.85),
the analyses have very high Q1 values equal to 0.98 (Case 1) and 0.96 (Case 2), respectively.
Figure 8a compares the calculated NRMSE values over the considered cardiac cycle. For
all time points, the PC-MRI data has a noticeably higher NRMSE than both DA analyses.
When comparing the two different data assimilation experiments, the time-independent anal-
ysis flow field shows a better agreement with the PIV measurements. Although the RMSEs
have been normalized by the velocity, higher values are still observed at time points that are
associated to higher peak flows (see also Figure 4). The Q1 values calculated in Figure 8b
support the previous discussion. The Q1 values are the lowest (i.e., poorest) for the raw
PC-MRI data, while they are almost always maximal (i.e., best) for the time-independent
urn al P
Due to recent advances in computational resources as well as medical imaging techniques,
the patient-specific investigation of intracranial blood flows has become increasingly popu-
lar [30, 31]. Limitations in numerical and experimental methods impair the acceptance of
intracranial flow investigations by the clinical community [8, 9, 38]. By introducing the tran-
sient data assimilation approach a better prediction of the velocity field becomes possible; this
has been validated by comparison with high-resolution PIV data. Previous studies have al-
ready tried to assimilate PC-MRI data and fluid dynamics simulations [20, 21, 22, 23, 26, 27];
however, validation was difficult due to the lack of a ground truth – or sometimes of any
real measurement data. A qualitative comparison of the velocity magnitude reveals a high
similarity between the ground-truth PIV data and the calculated DA analyses. In both
cases, the noise is reduced in comparison to the raw MRI data. The time-independent data
assimilation approach results in smooth velocity flow fields and continuous ensemble trajec-
tories throughout the whole cardiac cycle. From the calculated analysis flow field, the WSS
distribution can be calculated successfully, as shown in Figure 5. These observations are
further supported by the quantitative results. The spread of the velocity values in Figure 6
around the optimal red line is reduced from σ = 0.014 to σ = 0.01. At higher velocity values
a systematic error appears in the PC-MRI measurement data, which increases the velocity
differences of the corresponding data points. The LETKF can only account for random-
based error sources, but cannot improve systematic deviations. Therefore, the calculated
DA analyses are still biased by the systematic error sources of the original measurement
data. Nevertheless, NRMSE and Q1 values indicate that the analysis flow fields are much
closer to the ground truth than the raw PC-MRI measurement data. The noise reduction in-
dicates a better agreement with the general conservation laws (mass and momentum), while
PC-MRI data do not have any physical constraints implemented in the acquisition sequence.
In previous studies, the error sources for 4D flow data have been mainly attributed to
temporal and spatial downsampling, as well as acquisition artifacts such as partial volume
effects . For our own PC-MRI measurement data, we hypothesize a relative dependency
of the errors from the absolute velocity, which could perhaps be reduced by introducing
a scaling factor. We further assume that the value of this factor might be associated to
gradient inhomogeneities or to the Lorentz force. These effects possibly influence the phase
of the acquired MR signal and thus the quantitative flow data. Currently, a follow-up study
is being performed to quantify such possible artifacts and to identify the physical origin of
such a scaling factor.
urn al P
The spatial resolution can be increased for both, time-dependent and time-independent
data assimilation. However, the need and possibility for a temporal improvement remains
unclear. To calculate the pressure and WSS field, an increase in temporal resolution is not
really needed, but an improved spatial resolution is crucial. In terms of rupture or growth
prediction, it could be sufficient to only improve the spatial resolution and physical correct-
ness at selected time-points. Thus, the time-independent data assimilation approach should
currently be preferred, as possible parallel analyses calculations further reduce the compu-
tational time. Additionally, continuous ensemble inflow trajectories throughout the whole
cardiac cycle for all ensemble members lead to a smoother ensemble solution and smooth
analysis flow fields. With the introduction of a simplified cardiac model, as it was attempted
in the time-dependent approach, abrupt changes in the inflow trajectory after an analysis
calculation are introduced. This might distort the forward propagation of the individual en-
semble members. To ensure a proper reconstruction of the cardiac cycle trajectory, a better
cardiac model forecast must be introduced.
The computational time needed to calculate one CFD simulation is ≈ 32 CPU hours. This
is about the same time needed to calculate one analysis time step, without any paralleliza-
tion. By distributing the calculation of local analyses on 6 nodes, the overall computational
time was reduced by a factor of 4. The localization step not only enables parallel compu-
tations, but also reduces the ensemble size from 50 to 25 members in comparison to that
needed in . Fine-tuning of parameters such as the localization radius or ensemble size
will be considered in further studies. After obtaining an assimilated high-resolution velocity
field, clinically relevant parameters like pressure distribution and WSS can be derived.
urn al P
Limitations of the Study. The whole process aims to point out the feasibility of the
LETKF to merge CFD and PC-MRI data together. Therefore, several assumptions have
been made to reduce the complexity of the study, mainly with respect to the numerical
model and measurement data. Fluid-structure interactions are not considered, because the
influence of the time-dependent pressure on the walls of the silicone phantom appears visually
to be negligible. However, it might be important to consider the detection of the vessel
walls in future studies, in particular when using a more flexible material or direct phantom
models with thin walls. Patient-specific configurations have not been considered in this
study. However, the use of an idealized aneurysm model ensures well-controlled laminar
flow and, therefore, a better validation of the data assimilation step. The complexity of the
data assimilation steps has been further reduced by using underlying PC-MRI measurement
data that have a relatively high resolution, in order to achieve reliable state estimates for
the data assimilation step. Future studies will systematically decrease the resolution of the
measurement data to reduce scan times, but at the same time ideally keep reliable state
estimates with the help of data assimilation. Finally, as only a single configuration has been
considered in this initial study, it is obvious that the present findings should be confirmed
by later systematic investigations involving a sufficient number of patient-specific vessel
The current study describes an implementation of LETKF to merge PC-MRI measure-
ment data together with numerical flow simulations of an artificial intracranial aneurysm
geometry. A unique validation is ensured by high-resolution PIV measurements of the same
geometry. Both qualitative and quantitative validation demonstrate that the data assimila-
tion step works as a filter and increases the spatial resolution while reducing the noise of the
measurements. Conservation laws are enforced by DA onto the measurement data by the
Navier-Stokes equations, which increases the physical correctness of the data assimilation
results. Using a simplified cardiac forecast model, a forecast step in the classical sense has
been introduced to reconstruct the inflow conditions for the ensemble propagation. However,
the underlying PC-MRI measurements apparently show a systematic deviation of the true
values for high velocities; this systematic bias cannot be improved by LETKF. In order to
further improve the data assimilation results, a bias correction of the measurement data is
crucial. These findings emphasize the need for an improved quantification of uncertainty
for the underlying PC-MRI measurements. In this way, the observation error covariance
matrix can also be adjusted to better account for uncertainties in the measurement data.
Additionally, an improved model for the cardiac cycle could enhance the forecast of the
aneurysm flow in-between available measurement time steps. Finally, uncertain geometric
boundaries should be additionally taken into account within the data assimilation step for
future applications. Including fluid-structure interaction in the background model for future
patient-specific simulations would be a straightforward implementation, as this procedure
is already established in hemodynamic CFD simulations. Furthermore, PC-MRI is able to
urn al P
detect wall movements, even though the corresponding uncertainty will probably be in the or-
der of the spatial resolution. Accordingly, the determination of the true geometric boundary
after the assimilation step is one of the main tasks for the use of LETKF in patient-specific
applications. One possibility would consist of detecting the zero-contour of the calculated
analysis, considering that the near-wall velocities need to be zero.
To increase the reliability of the data assimilation algorithm a follow-up study will con-
sider systematic changes of the spatial and temporal resolution of the underlying measure-
ment data. The goal is to enable lower resolution PC-MRI measurements to reduce scan
times, while keeping at the same time reliable state estimates thanks to data assimilation.
Fine tuning of current data assimilation parameters is also important.
This work was in part conducted within the context of the International Graduate
School MEMoRIAL at Otto von Guericke University (OVGU) Magdeburg, Germany, kindly
supported by the European Structural and Investment Funds (ESF) under the program
“Sachsen-Anhalt WISSENSCHAFT Internationalisierung” (project no. ZS/2016/08/80646).
Fruitful discussions with the research team of Professor Sebastian Reich (University of Pots-
dam), in particular with Dr. Sahani Pathiraja, is warmly acknowledged.
Conflict of interest
The authors declare no conflict of interest. References
urn al P
 J. R. Cebral, M. Vazquez, D. M. Sforza, G. Houzeaux, S. Tateshima, E. Scrivano,
C. Bleise, P. Lylyk, C. M. Putman, Analysis of hemodynamics and wall mechan-
ics at sites of cerebral aneurysm rupture, J Neurointerv Surg 7 (7) (2015) 530–536.
 G. Janiga, P. Berg, S. Sugiyama, K. Kono, D. A. Steinman, The Computational Fluid
Dynamics Rupture Challenge 2013 – Phase I: Prediction of rupture status in intracranial
aneurysms, AJNR Am. J. Neuroradiol. 36 (3) (2015) 530–536. doi:10.3174/ajnr.A4157.
 J. Liu, J. Xiang, Y. Zhang, Y. Wang, H. Li, H. Meng, X. Yang, Morphologic and hemo-
dynamic analysis of paraclinoid aneurysms: Ruptured versus unruptured, J Neurointerv
Surg 6 (9) (2014) 658–663. doi:10.1136/neurintsurg-2013-010946.
 I. G. H. Jansen, J. J. Schneiders, W. V. Potters, P. van Ooij, R. van den Berg,
E. van Bavel, H. A. Marquering, C. B. L. M. Majoie, Generalized versus patient-
specific inflow boundary conditions in computational fluid dynamics simulations of cere-
bral aneurysmal hemodynamics, AJNR Am. J. Neuroradiol. 35 (8) (2014) 1543–1548.
 C. Karmonik, C. Yen, O. Diaz, R. Klucznik, R. G. Grossman, G. Benndorf, Tempo-
ral variations of wall shear stress parameters in intracranial aneurysms–importance of
patient-specific inflow waveforms for CFD calculations, Acta Neurochir 152 (8) (2010)
urn al P
 J. Jiang, C. Strother, Computational fluid dynamics simulations of intracranial
aneurysms at varying heart rates: A patient-specific study, J. Biomech. Eng. 131 (9)
(2009) 091001 (11 pages). doi:10.1115/1.3127251.  D. F. Kallmes, Point: CFD–computational fluid dynamics or confounding factor dis-
semination, AJNR Am. J. Neuroradiol. 33 (3) (2012) 395–396. doi:10.3174/ajnr.A2993.
 J. R. Cebral, H. Meng, Counterpoint: Realizing the clinical utility of Computational
Fluid Dynamics – Closing the gap, AJNR Am. J. Neuroradiol. 33 (3) (2012) 396–398.
 C. M. Strother, J. Jiang, Intracranial aneurysms, cancer, x-rays, and computational fluid dynamics, AJNR Am. J. Neuroradiol. 33 (6) (2012) 991–992. doi:10.3174/ajnr.A3163.
 M. Markl, A. Frydrychowicz, S. Kozerke, M. Hope, O. Wieben, 4D flow MRI, J. Magn. Reson. Imaging 36 (5) (2012) 1015–1036. doi:10.1002/jmri.23632.  J. Lotz, C. Meier, A. Leppert, M. Galanski, Cardiovascular flow measurement with
phase-contrast MR imaging: Basic facts and implementation, Radiographics: A review
publication of the Radiological Society of North America, Inc 22 (3) (2002) 651–671.
 N. Pelc, R. Herfkens, A. Shimakawa, Phase contrast cine magnetic resonance imaging, Magnetic Resonance Quarterly 7 (1991) 229–254.
 C. Roloff, D. Stucht, O. Beuing, P. Berg, Comparison of intracranial aneurysm
flow quantification techniques: Standard PIV vs stereoscopic PIV vs tomographic
PIV vs phase-contrast MRI vs CFD, J Neurointerv Surg 11 (3) (2019) 275–282.
urn al P
 J. Busch, D. Giese, L. Wissmann, S. Kozerk, Reconstruction of divergence-free velocity
fields from cine 3D phase-contrast flow measurements, Magn Reson Med 69 (2012)
 M. F. Sereno, B. K¨ohler, B. Preim, Comparison of divergence-free filters for cardiac 4D
PC-MRI data, in: Bildverarbeitung f¨ ur die Medizin 2018, Springer Berlin Heidelberg,
2018, pp. 139–144.
 P. Berg, D. Stucht, G. Janiga, O. Beuing, O. Speck, D. Th´evenin, Cerebral blood flow in
a healthy Circle of Willis and two intracranial aneurysms: Computational fluid dynamics
versus four-dimensional phase-contrast magnetic resonance imaging, J. Biomech. Eng.
136 (4) (2014) 041003 (9 pages). doi:10.1115/1.4026108.
 M. Lemke, J. Sesterhenn, Adjoint-based pressure determination from PIV data in com-
pressible flows – Validation and assessment based on synthetic data, Eur. J. Mech.
B/Fluids. 58 (2016) 29–38. doi:10.1016/j.euromechflu.2016.03.006.
 M. Lemke, J. Reiss, J. Sesterhenn, Pressure estimation from PIV like data of compress-
ible flows by boundary driven adjoint data assimilation, in: AIP Conference Proceed-
ings, 2016. doi:10.1063/1.4951773.  M.
 M. D’Elia, L. Mirabella, T. Passerini, M. Perego, M. Piccinelli, C. Vergara, A. Veneziani,
Applications of variational data assimilation in computational hemodynamics, in:
D. Ambrosi, A. Quarteroni, G. Rozza (Eds.), Modeling of Physiological Flows, Springer
Milan, 2012, pp. 363–394.
 M. D’Elia, M. Perego, A. Veneziani, A variational data assimilation procedure for the
incompressible Navier-Stokes equations in hemodynamics, J. Sci. Comput. 52 (2) (2012)
urn al P
 M. D’Elia, A. Veneziani, Uncertainty quantification for data assimilation in a steady in-
compressible Navier-Stokes problem, ESAIM Math. Model. Numer. Anal. 47 (4) (2013)
 S. W. Funke, M. Nordaas, Ø. Evju, M. S. Alnæs, K. A. Mardal, Variational data
assimilation for transient blood flow simulations: Cerebral aneurysms as an illustrative
example, Int J Numer Meth Bio. 35 (1) (2018) e3152.  X. Gao, Y. Wang, N. Overton, M. Zupanski, X. Tu, Data-assimilated computational
fluid dynamics modeling of convection-diffusion-reaction problems, J. Comput. Sci. 21
(2017) 38–59. doi:10.1016/j.jocs.2017.05.014.
 M. C. Rochoux, B. Cu´enot, S. Ricci, A. Trouv´e, B. Delmotte, S. Massart, R. Paoli,
R. Paugam, Data assimilation applied to combustion, Comptes Rendus M´ecanique
341 (1) (2013) 266–276. doi:10.1016/j.crme.2012.10.011.
 A. Bakhshinejad, V. L. Rayz, R. M. D’Souza, Reconstructing blood velocity profiles
from noisy 4D-PCMR data using ensemble Kalman filtering, in: Biomedical Engineering
Society (BMES) - Annual Meeting, Minneapolis, Minnesota, 2016.  A. Bakhshinejad, A. Baghaie, A. Vali, D. Saloner, V. L. Rayz, R. M. D’Souza,
Merging computational fluid dynamics and 4D Flow MRI using proper or-
thogonal decomposition and ridge regression, J. Biomech. 58 (2017) 162–173.
 G. Janiga, Quantitative assessment of 4D hemodynamics in cerebral aneurysms using
 G. Evensen, Data Assimilation - The Ensemble Kalman Filter, 2nd Edition, Springer International Publishing, 2009.
 M. Markl, F. Chan, M. Alley, K. Wedding, M. Draney, C. Elkins, D. Parker, R. Wicker,
C. A. Taylor, R. J. Herfkens, Time-resolved three-dimensional phase-contrast MRI, J.
Magn. Reson. Imaging 17 (4) (2003) 499–506. doi:10.1002/jmri.10272.
urn al P
 M. Markl, A. Harloff, T. Bley, M. Zaitsev, B. Jung, E. Weigang, M. Langer, J. Hennig,
A. Frydrychowicz, Time-resolved 3D MR velocity mapping at 3T: Improved navigator-
gated assessment of vascular anatomy and blood flow, J. Magn. Reson. Imaging 25 (4)
(2007) 824–831. doi:doi.org/10.1002/jmri.20871.  J. Bock, B. Kreher, J. Hennig, M. Markl, Optimized pre-processing of time-resolved
2D and 3D phase contrast MRI data, in: Proceedings of the 15th Annual Meeting of
ISMRM, Berlin, Germany, 2007, p. 3135.
 P. Berg, S. Saalfeld, S. Voß, O. Beuing, G. Janiga, A review on the reliabil-
ity of hemodynamic modeling in intracranial aneurysms: Why computational fluid
dynamics alone cannot solve the equation, Neurosurg. Focus 47 (1) (2019) E15.
 B. R. Hunt, E. J. Kostelich, I. Szunyogh, Efficient data assimilation for spatiotemporal
chaos: A local ensemble transform Kalman filter, Physica D: Nonlinear Phenomena
230 (1) (2007) 112–126. doi:10.1016/j.physd.2006.11.008.  J. Harlim, B. R. Hunt, Four-dimensional local ensemble transform Kalman filter: nu-
merical experiments with a global circulation model, Tellus A: Dynamic Meteorology
and Oceanography 59 (5) (2007) 731–748. doi:10.1111/j.1600-0870.2007.00255.x.
 E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E. Kalnay,
D. J. Patil, J. A. Yorke, A local ensemble Kalman filter for atmospheric data as-
similation, Tellus A: Dynamic Meteorology and Oceanography 56 (5) (2004) 415–428.
 C. H. Bishop, B. J. Etherton, S. J. Majumdar, Adaptive sampling with the ensemble
transform Kalman filter. Part I: Theoretical aspects, Mon. Weather Rev. 129 (3) (2001)
urn al P
 P. Berg, C. Roloff, O. Beuing, S. Voss, S.-I. Sugiyama, N. Aristokleous, A. S. Anayiotos,
N. Ashton, A. Revell, N. W. Bressloff, A. G. Brown, B. Jae Chung, J. R. Cebral,
G. Copelli, W. Fu, A. Qiao, A. J. Geers, S. Hodis, D. Dragomir-Daescu, E. Nordahl,
Y. Bora Suzen, M. Owais Khan, K. Valen-Sendstad, K. Kono, P. G. Menon, P. G. Albal,
O. Mierka, R. M¨ unster, H. G. Morales, O. Bonnefous, J. Osman, L. Goubergrits, J. Pal-
lares, S. Cito, A. Passalacqua, S. Piskin, K. Pekkan, S. Ramalho, N. Marques, S. Sanchi,
ˇ K. R. Schumacher, J. Sturgeon, H. Svihlov´ a, J. Hron, G. Usera, M. Mendina, J. Xiang,
H. Meng, D. A. Steinman, G. Janiga, The Computational Fluid Dynamics Rupture
Challenge 2013 – Phase II: Variability of hemodynamic simulations in two intracranial
aneurysms, J. Biomech. Eng. 137 (121008 (13 pages)). doi:10.1115/1.4031794.
 P. van Ooij, A. Gu´edon, C. Poelma, J. Schneiders, M. C. M. Rutten, H. A. Marquer-
ing, C. B. Majoie, E. vanBavel, A. J. Nederveen, Complex flow patterns in a real-size
intracranial aneurysm phantom: phase contrast mri compared with particle image ve-
locimetry and computational fluid dynamics, NMR in Biomedicine 25 (1) (2012) 14–26.
urn al P