# Exploring Probabilistic Distance Fields in Robotics

Lan Wu. Email<EMAIL_ADDRESS>
Robotics Institute, University of Technology Sydney

##  I Introduction

A Robot’s successful performance of a highly intelligent mission typically
entails the integration of several foundational research tasks. Each task
imposes unique requirements on the task-oriented representation, reflecting
robotic representations’ diversity and complexity. For example, feature-based
representations [1, 2, 3] are widely used by most SLAM methods for
localisation and mapping tasks [4, 5] and by computer vision for the context
of Structure from Motion [6, 7]. However, it is difficult to reconstruct a
dense map from sparse features for visualisation tasks and to use sparse
landmarks for planning and human interaction. Raw and dense representations
provide decent visual information and high texture in geometry via a large
point cloud [8, 9]. The drawback is they usually require managing and storing
a large amount of data. Moreover, they give a low-level understanding of the
scene, which is of limited applicability for navigation tasks. Volumetric
Dense Representations [10, 11, 12, 13, 14, 15] such as Truncated Signed
Distance Field (TSDF) are more suitable for planning, obstacle avoidance,
manipulation, and reconstruction. However, without matching sparse features to
identify correspondences, it is difficult to perform a consistent alignment
for localisation [16].

Thus, designing a representation for every different task of a complex mission
is impractical and costly. Unified and task-centric representations that
fulfil multiple requirements remain unexplored, demanding more innovation to
shape the future of robotic intelligence. The most relevant representations
that aim to cater for multiple tasks are as follows. Compared to TSDF,
Euclidean Signed Distance Field (ESDF) [17, 18, 19, 20] contains distance even
with gradient fields over the entire exploring environment suitable for
localisation, mapping, planning and manipulation. A major challenge is that it
is expensive to compute and propagate the distance information over the whole
space when a task might only require this information at discrete sparse
points. Recent developments with neural networks include a continuous implicit
representation for reconstruction [21, 22, 23, 24, 23, 25]. Such neural
representations have shown great promise for mapping, localisation [26, 27,
28] as well as planning [29, 30, 31, 32]. Although these neural
representations serve as a sound representation for multiple tasks, they are a
‘black-box’ which requires further theoretical formulation.

On the other hand, Gaussian Process (GP) [33, 34, 35, 36, 37] probabilistic
representation is a well-established approach that deals with the sparseness
of sensory information, manages and handles uncertainty, allows a continuous
prediction and even supports high-level semantic understanding of the scene. I
believe GP-based probabilistic representations have the potential to offer a
theoretical sound, universal and ”one for all” representation solution for
localisation, visualisation, mapping, motion planning, manipulation, human-
robot interaction and many more demands in robotics.

As shown in Fig. 1, my work proposes a GP-based probabilistic distance field
(GPDF) representation that mathematically models the fundamental property of
Euclidean distance field (EDF) along with gradients, surface normals and dense
reconstruction [38]. EDF with gradients is crucial for motion planning and
dynamic obstacle avoidance in human-robot interaction tasks. Surface normal
directly feeds tasks such as grasping. A densely reconstructed mesh can be
used for visualisation and graphical interfaces. The accurately inferred
distance fields allow pose estimation via direct optimisation of the
alignment. The distance and gradient are only computed when tasks are querying
the representation. Furthermore, the nature of the GP representation provides
uncertainty for probabilistic data fusion. Overall, GPDF overcomes the
limitations of other state-of-the-art representations in robotics, as we
highlighted above. GPDF offers a unified solution that transcends the need to
design multiple representations for diverse underlying tasks.

I anticipate that GPDF will form the basis representation of the full robot
autonomy. GPDF can be extended for efficient online running, more robust
performance, for tightly integrating with a new motion planner taking
advantage of all the information output from GPDF, and for applying it to the
human-robot collaboration and manipulation tasks. GPDF focus on spatial-only
information at the current stage. Future vision considers expanding it to
include temporal information or even equipped with scene understanding. The
nature of the GP has the advantage of integrating GPDF with Gaussian spatting
for high-quality reconstruction. Currently, GPDF works with depth sensors. The
probabilistic nature of GPDF can be directly used for data fusion, combining
multiple sensors in the future for spatial, temporal and semantic information.
In what follows, I will present the progress to support the prospection
throughout my PhD and first-year postdoc.

Figure 1: Block Diagram of the proposed GP-based representation architecture
for multiple robotics applications and future research objectives.

##  II Progress to Date

###  II-A Probabilistic Distance Fields

I initially proposed GPDF as a novel representation [39]. My key contribution
is the realisation that the regularised Eikonal equation can be simply solved
by applying the logarithmic transformation to a GP formulation to approximate
the accurate Euclidean distance fields. The proposed approach does not require
post-processing or wave propagation to recover the distance fields. Our
experiments show that it produces the most accurate results for distance and
surface reconstruction. However, it requires a trade-off between accuracy and
interpolation abilities. These are limiting factors for specific applications
that require highly accurate values of the Euclidean distance fields, such as
in the case of echolocation for robotic inspection. Additionally, due to the
non-linear transformation that is applied over the GP fields, the variance
prediction of our proposed method makes it hard to model the field’s
uncertainty reliably. Following the theoretical formulation in [39], we have
introduced the so-called reverting GP distance field [40] that increased the
accuracy and provided a proxy for its uncertainty. GP-based distance fields
are appealing as they can represent complex environments with non-parametric
models. The accuracy of the proposed distance fields enables us to tackle the
multiple key problems in robotics with noisy measurements and uncertainty.

###  II-B Unified Representation

My work proposed GPDF as a unified representation for the multiple robotic
tasks in [41]. It allows us to use a single representation for localising,
mapping, and path planning. The continuous and probabilistic nature provides
an accurate dense surface for mapping. The accurately inferred distance fields
allow pose estimation via direct optimisation of the alignment between the
global map representation and a new set of measurements. Furthermore, GPDF is
integrated with a optimisation-based planners.

###  II-C Towards Efficiency of the Probabilistic Representation

I have been working on efficient GP modelling. Firstly, I proposed a novel
submapping method [42], following the topology of the scene to generate
conditional independent submaps of GP. Large datasets can efficiently be
processed by the proposed pipeline, producing not only a surface but also the
uncertainty information of the representation. Secondly, the pseudo inputs
applied in efficient GP frameworks are a subset of the full data in general.
We proposed an alternative solution [43] to constrain the locations of pseudo
data using the distance field itself provided by our GP kernel functions. We
demonstrate the effectiveness and sound performance of our method in different
applications.

##  III Ongoing and Future Work

Developing and maintaining multiple representations tailored to each task
requires significant time, computational resources, and expertise. In
contrast, adopting a unified GPDF representation approach offers a more
streamlined and cost-effective solution. I believe that GPDF serves as the
cornerstone for robots to accomplish more complex and challenging tasks. By
leveraging GPDF, robots can navigate through intricate environments,
understand spatial relationships, and interact with objects and humans
seamlessly. Toward this objective, some ongoing and future works are
highlighted.

###  III-A Online Operation of Probabilistic Distance Fields

The main concern with GP-based algorithms is the cubic computational
complexity [33]. Substantial efforts have been made in the literature and my
progress to date. However, a notable gap remains to achieve the desired online
performance. I am developing an online GPDF framework that provides distance,
gradients, repulsive information, rapid continuous reconstruction and colour
as output. Our framework develops a fast access data structure to maintain
feasible, long-term and large-sized mapping and exploration. It is important
to have online GPDF since most robotic tasks have to be conducted in real-
time. Our ongoing progress (in Fig. 2a) has successfully generated varying
information in real-time and is more efficient than similar state-of-the-art
frameworks [17, 18].

###  III-B Capability for Human-Robot Collaboration

We are extending the GPDF for Human-Robot Collaboration (HRC) and grasping
tasks. One essential requirement for the perception of the interactive scene
is to deal with dynamic objects and human interaction in the area. My previous
research outcomes have yielded progress in dynamic application [44, 43].
However, further work is needed to bridge the remaining distance to attain the
interactive distance fields with effective dynamic capability and efficiency.
As shown in Fig. 2b, we have open-sourced our framework [45] for
HRC.333https://uts-ri.github.io/IDMP.

###  III-C Robust Performance for Localisation

I am developing the machine perception part for a project delivering a form of
technologically enhanced human echolocation that enables blind users to
perceive their surroundings in precise detail. One requirement is to solve the
well-known SLAM problem [46, 47]. One of my future works is to combine
feature-based odometry with our GPDF-based odometry [41] to create a more
robust localisation framework. The solution needs to be computationally light
and robust, enabling it to run smoothly in an embedded device with AI glasses.

(a)

(b)

Figure 2: a) Online GPDF framework performance for Sec. III-A. b) Dynamical
GPDF for Human-robot collaboration application for Sec. III-B.

###  III-D Theoretical New Planner Based on GP Distance Fields

With various information directly queried from the GPDF, other applications
will take advantage of the novel representation solution. We are working on
designing a new planning algorithm with distance, gradient and repulsive
information.

###  III-E Vision for Future Work

Currently, the primary focus of GPDF lies in spatial representation. However,
our future vision entails expanding its capabilities to include temporal
information and scene understanding. The nature of GP has the potential to
integrate with Gaussian Splatting for radiance field rendering and be used for
probabilistic data fusion in future applications.

UTS

    University of Technology Sydney
RI

    Robotics Institute
FEIT

    Faculty of Engineering and Information Technology
CAS

    Centre for Autonomous Systems
WiEIT

    Women in Engineering and IT
1D

    One-Dimensional
2D

    Two-Dimensional
2.5D

    Two-and-a-Half-Dimensional
3D

    Three-Dimensional
GP

    Gaussian Process
GPIS

    Gaussian Process Implicit Surfaces
MAP

    Maximum A Posteriori
MLE

    Maximum Likelihood Estimation
ICP

    Iterative Closest Point
MVE

    Multi-View Environment
OG

    Occupancy Grid
CHOMP

    Covariant Hamiltonian Optimization for Motion Planning
FIESTA

    Fast Incremental Euclidean DiSTAnce Fields
GPUs

    Graphics Processing Units
KD

    K-Dimensional
Log-GPIS

    Log-Gaussian Process Implicit Surfaces
LiDAR

    Light Detection And Ranging Sensor
SLAM

    Simultaneous Localisation and Mapping
TDF

    Truncated Distance Field
EDF

    Euclidean Distance Field
SDF

    Signed Distance Field
TSDF

    Truncated Signed Distance Field
ESDF

    Euclidean Signed Distance Field
GPIS-SDF

    GPIS with signed distance function
RGB

    Red-Green-Blue
RGB-D

    Red-Green-Blue-Depth
RMSE

    Root Mean Squared Error
CI

    Conditional Independent
FP

    Forward Propagation
BP

    Backward Propagation
D-SKI

    Structured Kernel Interpolation framework with Derivatives
SKI

    Structured Kernel Interpolation method
D-SKI-CI-Fusion

    Structured Kernel Interpolation with Derivatives and Conditional Independent Fusion
MVMs

    Matrix-Vector Multiplications method
PDE

    Partial Differential Equation
Log-GPIS-MOP

    Log-Gaussian Process Implicit Surface for Mapping, Odometry and Planning
RANSAC

    Random Sample Consensus
Dynamic-GPDF

    Dynamic Gaussian Process Distance Fields

## References

  * [1] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” in _2011 International conference on computer vision_. Ieee, 2011, pp. 2564–2571.
  * [2] H. Bay, T. Tuytelaars, and L. V. Gool, “Surf: Speeded up robust features,” in _European conference on computer vision_. Springer, 2006, pp. 404–417.
  * [3] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International journal of computer vision_ , vol. 60, no. 2, pp. 91–110, 2004.
  * [4] R. Mur-Artal and J. D. Tardós, “ORB-SLAM2: an open-source SLAM system for monocular, stereo, and RGB-D cameras,” _IEEE Transactions on Robotics (T-RO)_ , vol. 33, no. 5, pp. 1255–1262, 2017.
  * [5] T. Qin, P. Li, and S. Shen, “Vins-mono: A robust and versatile monocular visual-inertial state estimator,” _IEEE Transactions on Robotics_ , 2018.
  * [6] O. Özyeşil, V. Voroninski, R. Basri, and A. Singer, “A survey of structure from motion*.” _Acta Numerica_ , vol. 26, pp. 305–364, 2017.
  * [7] S. Fuhrmann, F. Langguth, and M. Goesele, “Mve: A multi-view reconstruction environment,” in _Eurographics Workshop on Graphics and Cultural Heritage (GCH)_ , 2014, pp. 11–18.
  * [8] R. B. Rusu, Z. C. Marton, N. Blodow, M. Dolha, and M. Beetz, “Towards 3d point cloud based object maps for household environments,” _Robotics and Autonomous Systems_ , vol. 56, no. 11, pp. 927–941, 2008.
  * [9] M. Labbe and F. Michaud, “Online global loop closure detection for large-scale multi-session graph-based slam,” in _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2014, pp. 2661–2666.
  * [10] Y. Roth-Tabak and R. Jain, “Building an environment model using depth information,” _Computer_ , vol. 22, no. 6, pp. 85–90, 1989.
  * [11] H. Moravec, “Robot spatial perception by stereoscopic vision and 3d evidence grids,” _Perception_ , vol. 474, p. 475, 1996.
  * [12] B. Curless and M. Levoy, “A volumetric method for building complex models from range images,” in _Proc. of the Conference on Computer Graphics and Interactive Techniques_ , 1996, p. 303–312.
  * [13] S. Izadi, D. Kim, O. Hilliges, D. Molyneaux, R. Newcombe, P. Kohli, J. Shotton, S. Hodges, D. Freeman, A. Davison _et al._ , “Kinectfusion: real-time 3d reconstruction and interaction using a moving depth camera,” in _Proc. of the ACM symposium_ , 2011, pp. 559–568.
  * [14] R. Wagner, U. Frese, and B. Bäuml, “3d modeling, distance and gradient computation for motion planning: A direct gpgpu approach,” in _2013 IEEE International Conference on Robotics and Automation_. IEEE, 2013, pp. 3586–3592.
  * [15] R. Wagner, U. Frese, and B. Bauml, “Real-time dense multi-scale workspace modeling on a humanoid robot,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013.
  * [16] V. Reijgwart, A. Millane, H. Oleynikova, R. Siegwart, C. Cadena, and J. Nieto, “Voxgraph: Globally consistent, volumetric mapping using signed distance function submaps,” _IEEE Robotics and Automation Letters_ , 2020.
  * [17] H. Oleynikova, Z. Taylor, M. Fehr, R. Siegwart, and J. Nieto, “Voxblox: Incremental 3d euclidean signed distance fields for on-board mav planning,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 1366–1373.
  * [18] Y. Pan, Y. Kompis, L. Bartolomei, R. Mascaro, C. Stachniss, and M. Chli, “Voxfield: Non-projective signed distance fields for online planning and 3d reconstruction,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2022)_ , 2022.
  * [19] D. Zhu, C. Wang, W. Wang, R. Garg, S. Scherer, and M. Q.-H. Meng, “Vdb-edt: An efficient euclidean distance transform algorithm based on vdb data structure,” _arXiv preprint arXiv:2105.04419_ , 2021.
  * [20] L. Han, F. Gao, B. Zhou, and S. Shen, “Fiesta: Fast incremental euclidean distance fields for online motion planning of aerial robots,” in _2019 IEEE/RSJ IROS_. IEEE, 2019, pp. 4423–4430.
  * [21] A. Gropp, L. Yariv, N. Haim, M. Atzmon, and Y. Lipman, “Implicit geometric regularization for learning shapes,” _arXiv preprint arXiv:2002.10099_ , 2020.
  * [22] J. J. Park, P. Florence, J. Straub, R. Newcombe, and S. Lovegrove, “Deepsdf: Learning continuous signed distance functions for shape representation,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 165–174.
  * [23] B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, and R. Ng, “Nerf: Representing scenes as neural radiance fields for view synthesis,” _Communications of the ACM_ , pp. 99–106, 2021.
  * [24] V. Sitzmann, J. Martel, A. Bergman, D. Lindell, and G. Wetzstein, “Implicit neural representations with periodic activation functions,” _NeurIPS_ , vol. 33, pp. 7462–7473, 2020.
  * [25] L. Mescheder, M. Oechsle, M. Niemeyer, S. Nowozin, and A. Geiger, “Occupancy networks: Learning 3d reconstruction in function space,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 4460–4470.
  * [26] E. Sucar, S. Liu, J. Ortiz, and A. J. Davison, “imap: Implicit mapping and positioning in real-time,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 6229–6238.
  * [27] A. Rosinol, J. J. Leonard, and L. Carlone, “Nerf-slam: Real-time dense monocular slam with neural radiance fields,” _arXiv preprint arXiv:2210.13641_ , 2022.
  * [28] Z. Zhu, S. Peng, V. Larsson, W. Xu, H. Bao, Z. Cui, M. R. Oswald, and M. Pollefeys, “Nice-slam: Neural implicit scalable encoding for slam,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 12 786–12 796.
  * [29] J. Ortiz, A. Clegg, J. Dong, E. Sucar, D. Novotny, M. Zollhoefer, and M. Mukadam, “isdf: Real-time neural signed distance fields for robot perception,” _arXiv preprint arXiv:2204.02296_ , 2022.
  * [30] D. Driess, J.-S. Ha, M. Toussaint, and R. Tedrake, “Learning models as functionals of signed-distance fields for manipulation planning,” in _Conference on Robot Learning_. PMLR, 2022, pp. 245–255.
  * [31] M. Adamkiewicz, T. Chen, A. Caccavale, R. Gardner, P. Culbertson, J. Bohg, and M. Schwager, “Vision-only robot navigation in a neural radiance world,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4606–4613, 2022.
  * [32] M. Pantic, C. Cadena, R. Siegwart, and L. Ott, “Sampling-free obstacle gradients and reactive planning in neural radiance fields (nerf),” _arXiv preprint arXiv:2205.01389_ , 2022.
  * [33] C. E. Rasmussen and C. K. Williams, _Gaussian Processes for Machine Learning_. Cambridge, Mass.: MIT Press, 2006.
  * [34] J. Mahler, S. Patil, B. Kehoe, J. Van Den Berg, M. Ciocarlie, P. Abbeel, and K. Goldberg, “Gp-gpis-opt: Grasp planning with shape uncertainty using gaussian process implicit surfaces and sequential convex programming,” in _2015 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2015, pp. 4919–4926.
  * [35] W. Martens, Y. Poffet, P. R. Soria, R. Fitch, and S. Sukkarieh, “Geometric priors for gaussian process implicit surfaces,” _IEEE Robotics and Automation Letters (RA-L)_ , pp. 373–380, 2017.
  * [36] S. Dragiev, M. Toussaint, and M. Gienger, “Gaussian process implicit surfaces for shape estimation and grasping,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2011, pp. 2845–2850.
  * [37] S. Kim and J. Kim, “Building occupancy maps with a mixture of gaussian processes,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2012, pp. 4756–4761.
  * [38] L. Wu, “Probabilistic implicit surfaces for localisation, mapping and planning,” Ph.D. dissertation, University of Technology Sydney (Australia), 2023.
  * [39] L. Wu, K. M. B. Lee, L. Liu, and T. Vidal-Calleja, “Faithful euclidean distance field from log-gaussian process implicit surfaces,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2461–2468, 2021.
  * [40] C. Le Gentil, O.-L. Ouabi, L. Wu, C. Pradalier, and T. Vidal-Calleja, “Accurate gaussian-process-based distance fields with applications to echolocation and mapping,” _IEEE Robotics and Automation Letters_ , 2023.
  * [41] L. Wu, K. M. B. Lee, C. Le Gentil, and T. Vidal-Calleja, “Log-GPIS-MOP: A Unified Representation for Mapping, Odometry, and Planning,” _IEEE Transactions on Robotics_ , pp. 1–17, 2023. [Online]. Available: https://ieeexplore.ieee.org/document/10202666/
  * [42] L. Wu, R. Falque, V. Perez-Puchalt, L. Liu, N. Pietroni, and T. Vidal-Calleja, “Skeleton-based conditionally independent gaussian process implicit surfaces for fusion in sparse to dense 3d reconstruction,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1532–1539, 2020.
  * [43] L. Wu, C. Le Gentil, and T. Vidal-Calleja, “Pseudo inputs optimisation for efficient gaussian process distance fields,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2023, pp. 7249–7255.
  * [44] L. Liu, S. Fryc, L. Wu, T. L. Vu, G. Paul, and T. Vidal-Calleja, “Active and Interactive Mapping With Dynamic Gaussian Process Implicit Surfaces for Mobile Manipulators,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3679–3686, Apr. 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9361306/
  * [45] U. Ali, L. Wu, A. Mueller, F. Sukkar, T. Kaupp, and T. Vidal-Calleja, “Interactive distance field mapping and planning to enable human-robot collaboration,” _arXiv preprint arXiv:2403.09988_ , 2024.
  * [46] J. Aulinas, Y. Petillot, J. Salvi, and X. Lladó, “The slam problem: a survey,” _Artificial Intelligence Research and Development_ , pp. 363–371, 2008.
  * [47] C. Cadena, L. Carlone, H. Carrillo, Y. Latif, D. Scaramuzza, J. Neira, I. Reid, and J. Leonard, “Past, present, and future of simultaneous localization and mapping: Towards the robust-perception age,” _Trans. on Rob._ , p. 1309–1332, 2016.

  *[TSDF]: Truncated Signed Distance Field
  *[ESDF]: Euclidean Signed Distance Field


# Physics-informed neural networks for blood flow inverse problems

Jeremías Garay  Jocelyn Dunstan  Sergio Uribe  Francisco Sahli Costabal
<EMAIL_ADDRESS>Department of Mechanical and Metallurgical Engineering,
Pontifical Catholic University of Chile  Center of Biomedical Imaging,
Pontifical Catholic University of Chile  Department of Medical Imaging and
Radiation Sciences, Monash University  Millennium Institute for Intelligent
Healthcare Engineering (iHealth), Chile  Department of Computer Science,
Pontifical Catholic University of Chile  Center for Mathematical Modeling,
University of Chile  Institute for Mathematical Computing Engineering,
Pontifical Catholic University of Chile  Institute for Biological and Medical
Engineering, Pontifical Catholic University of Chile  Department of Radiology,
Pontifical Catholic University of Chile

###### Abstract

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

###### keywords:

Physics-informed neural networks , hemodynamics , reduced-order modeling ,
blood flow , patient-specific model

††journal: Computer Methods in Applied Mechanics and Engineering

##  1 Introduction

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

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

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

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

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

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

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

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

##  2 Methods

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

###  2.1 The mathematical model

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

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

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

$\begin{cases}\displaystyle\rho\frac{\partial\textbf{u}}{\partial
t}+\rho\big{(}\textbf{u}\cdot\nabla\big{)}\textbf{u}-\mu\Delta\textbf{u}+\nabla
p=0\quad\text{in}\quad\Omega\times[0,T],\\\\[5.69046pt]
\nabla\cdot\textbf{u}=0\quad\text{in}\quad\Omega\times[0,T],\\\\[5.69046pt]
\textbf{u}=\textbf{u}_{inlet}(t)\quad\text{on}\quad\Gamma_{in}\times[0,T],\\\\[5.69046pt]
\textbf{u}=\textbf{0}\quad\text{on}\quad\Gamma_{wall}\times[0,T],\\\\[5.69046pt]
\displaystyle\mu\frac{\partial\textbf{u}}{\partial
n}-p\textbf{n}=-P_{k}(t)\displaystyle\textbf{n}\quad\text{on}\quad\Gamma_{k}\times[0,T],\quad
k=1,\dots,K\end{cases}$ (1)

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

$\begin{cases}P_{k}=R_{p,k}\ Q_{k}+\pi_{k},\\\\[5.69046pt]
Q_{k}=\displaystyle\int_{\Gamma_{k}}\textbf{u}\cdot\textbf{n},\\\\[5.69046pt]
\displaystyle
C_{k}\frac{d\pi_{k}}{dt}+\frac{\pi_{k}}{R_{d,k}}=Q_{k}.\end{cases}$ (2)

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

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

$\begin{cases}P_{k}=(R_{p,k}+R_{d,k})\ Q_{k}=R_{t,k}\ Q_{k},\\\\[5.69046pt]
Q_{k}=\displaystyle\int_{\Gamma_{k}}\textbf{u}\cdot\textbf{n}.\\\ \end{cases}$
(3)

with the distal and proximal resistances absorbed in a total resistance
parameter $R_{t}$.

Figure 1: Aortic domain $\Omega$. The Windkessel outlets are defined at
$\Gamma_{1}$ to $\Gamma_{5}$ consisting of 2 resistors ($R_{d}$, $R_{p}$) and
1 capacitor $C$ each.

###  2.2 Physics-Informed Neural Networks (PINNs)

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

$\big{(}\textbf{u}(\textbf{x},t),\
p(\textbf{x},t)\big{)}=\mathscr{NN}(\textbf{x},t;\mathbb{W},\textbf{b})$ (4)

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

$\mathscr{L}_{tot}=\mathscr{L}_{NS}+\mathscr{L}_{WK}+\mathscr{L}_{BC}+\mathscr{L}_{u,data}+\mathscr{L}_{p,data}+\mathscr{L}_{gradp}.$
(5)

where the individual loss components are defined as:

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

The last added term (11) is included to promote a null pressure gradient at
the Windkessel outlets since the coupled equations between the 3D and 0D
systems assumed the pressure to be constant at the interface. In all cases,
$||\cdot||$ represents the $L_{2}$ norm of the quantities. We make all
equations dimensionless for the loss function to obtain normalized predictions
of velocity and pressure, which improves the training procedure. We use the
following normalization of the physical quantities:

$\displaystyle\textbf{x}\longrightarrow\textbf{x}/L$
$\displaystyle\textbf{u}\longrightarrow\textbf{u}/U$ $\displaystyle
t\longrightarrow t/(L/U)$ $\displaystyle p\longrightarrow p/(\rho U^{2})$

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

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

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

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

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

Figure 2: Feed-forward neural network architecture used for the transient
problem. The network represents a vector potential $\Phi$ and the pressure
field. After physics-informed training, the output of the network is used for
reconstruct the 3D velocity field and the Windkessel parameters of all the
domain outlets.

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

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

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

  * i)

We use a vector potential representation of the velocity u to hard-impose mass
conservation of the solution. Consequently, the velocity is assumed to be
originated from the curl of a vector potential $\Phi(\textbf{x},t)$ [43] as:

$\textbf{u}(\textbf{x},t)=\nabla\times\Phi(\textbf{x},t).$ (13)

With this definition, the velocity field $vec{u}$ is always divergence free.
We find that this change of variable improves the quality of the results and
convergence of the method at the expense increasing the computational time.

  * ii)

We tune each weight of the total loss function automatically following the
algorithm presented in [44], using the $L_{1}$ norm of the gradient of the
physics loss. We perform the update every 10 epochs. The new values of the
weights are exponentially averaged with their respective history following the
rule:

$\lambda_{new,i}=(1-\alpha)\cdot\lambda_{old,i}+\alpha\cdot\frac{||\nabla\mathscr{L}_{phys}||_{1}}{\overline{\nabla\mathscr{L}_{i}}}\quad\quad\forall
i=1,2,...,N_{L}$

being $N_{L}$ the number of loss terms. In all cases, we set $\alpha$ to
$0.1$.

  * iii)

We divide the training into two stages. First, we randomly initialize the
Windkessel model parameters $R_{p,k},R_{d,k},C_{k}$ and fix them. Then, we
train only the network parameters $\mathbb{W},\textbf{b}$. In a second stage,
we allow the optimizer to also change the Windkensell model parameters. We
found that splitting the training in this way leads to better results since
the network had time first to adapt the representation of physical velocities
and pressures and then to estimate the Windkessel parameters.

  * iv)

We compute the compliance parameter $C$ from the distal resistance $R_{d}$ and
an estimation of the decay time $R_{d}\cdot C$. This time was obtained at the
beginning of the training from the mean pressure curve used as an extra
measurement. Furthermore, the decay time was computed by fitting an
exponential curve from the valve closure time and the end of the diastole, an
interval at which the Windkessel model assumes a decay behavior modulated by
$\sim\exp(-t/R_{d}C)$. The decay time was assumed constant for all outlets at
the aorta.

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

$\mathbb{W},\textbf{b}=\arg\min_{W,b}\ \mathscr{L}_{tot}(W,b).$ (14)

|  |  Steady Problem |  Transient Problem |
---|---|---|---|---
Estimation |  |  |  u , $p$ , $\big{\\{}R_{tot}\big{\\}}^{5}_{k=1}$ |  |  u , $p$ , $\big{\\{}R_{p},R_{d},C_{d}\big{\\}}^{5}_{k=1}$ |
Table 1: Different quantities to estimate in the steady and transient
problems.

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

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

Hyperparameter |  Steady Problem |  Transient Problem
---|---|---
Batchsize |  |  $1580$ |  | $3000$
Learning Rate |  |  $10^{-3}$ |  | $10^{-3}$
Learning Rate of Parameters |  |  $10^{-2}$ |  | $10^{-2}$
Learning Rate Scheduler |  | Yes |  | Yes
Hidden layers |  |  $7$ |  | $7$
No. neurons per layer |  |  $220$ |  | $220$
Hidden layers ($\pi$) |  |  $-$ |  | $6$
No. neurons per layer ($\pi$) |  |  $-$ |  | $10$
Activation Function |  | swish |  | swish
Epochs |  |  $1250$ |  | $120+1250$
Potential Vector Representation |  | No |  | Yes
Normalized Length ($L$) |  |  $1$ cm |  |  $0.5$ cm
Normalized Velocity ($U$) |  |  $5000$ cm/s |  |  $120$ cm/s
Loss weights |  | Manual |  | Automatic
| $\lambda_{phys}$ | 1.5 |  | $-$
| $\lambda_{data,u}$ | 1.0 |  | $-$
| $\lambda_{data,p}$ | 1.0 |  | $-$
| $\lambda_{BC}$ | 6.0 |  | $-$
| $\lambda_{windk}$ | 1.0 |  | $-$
Table 2: Hyperparameters used for training the PINNs

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

###  2.3 The Reference Solution

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

$\begin{cases}\displaystyle\textbf{u}=\textbf{0},\\\ \pi_{k}=78.8\
\rm{mmHg}\quad\text{for}\quad k=1,\dots,K.\end{cases}$

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

$\textbf{u}(\textbf{x},t)_{inlet}=\textbf{u}_{stokes}(\textbf{x})f(t)\
\textbf{n},$ (15)

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

$A_{inlet}=\gamma\int_{\Gamma_{inlet}}(\textbf{u}-\textbf{u}_{inlet})\cdot\textbf{v}$
(16)

with the parameter $\gamma$ fixed in $10^{5}\ \rm{gr/(cm^{2}\cdot s)}$.

(a)

(b)

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

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

Parameter | $\Gamma_{1}$ | $\Gamma_{2}$ | $\Gamma_{3}$ | $\Gamma_{4}$ | $\Gamma_{5}$
---|---|---|---|---|---
$R_{p}\ (\rm{dyn\cdot s\cdot cm^{-5}})$ | 713 | 713 | 602 | 689 | 98
$R_{d}\ (\rm{dyn\cdot s\cdot cm^{-5}})$ | 12023 | 12023 | 10143 | 11609 | 1650
$C_{d}\ (\rm{dyn^{-1}\cdot cm^{5}})$ | 8.256e-5 | 8.256e-5 | 9.785e-5 | 8.55e-5 | 6.015e-4
Table 3: Three-element Windkessel parameters for every outlet.

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

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

(a)

(b)

Figure 4: Reference simulation results. The blood flow velocity at peak
systole and its flow rates and average pressure curves are depicted in (a) and
(b), respectively.

###  2.4 PC-MRI measurement generation

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

$\displaystyle M^{u}_{meas}$ $\displaystyle=$ $\displaystyle M_{0}\
\exp(i\phi_{0}+i\pi u_{MRI}/venc),$ (17) $\displaystyle M^{0}_{meas}$
$\displaystyle=$ $\displaystyle M_{0}\ \exp(i\phi_{0}),$ (18)

where $M^{u}_{meas}$ is the velocity-encoded magnetization and $M^{0}_{meas}$
is an extra measurement usually done to capture the reference phase
$\phi_{0}$. The interpolated velocity is represented by $u_{MRI}$ while
$M_{0}$ is assumed constant within the vessel’s lumen. Moreover, Gaussian
noise is added to the magnetization components, producing a fixed signal-to-
noise ratio of $18$ dB. Consequently, a non-Gaussian noise distribution is
induced in the reconstructed velocity, as is generally observed in PC-MRI
measurements [50]. The final velocity field is then reconstructed using the
following equation:

$\textbf{u}_{meas}=venc\
\frac{\angle\big{(}M^{u}_{meas}/M^{0}_{meas}\big{)}}{\pi},$ (19)

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

(a)

(b)

Figure 5: Synthetic PC-MRI measurements obtained from the simulation. The
axial planes taken from the aorta are shown in (a), while the resulting PC-MRI
velocity measurements are shown in from frontal view in (b).

##  3 Results

###  3.1 Steady problem

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

Figure 6 shows the evolution of the estimated total resistance by every
training epoch and Windkessel outlet. We normalize all values to ease the
visualization of the parameter evolution. It can be seen that all parameters
stabilize rapidly to their reference value, yet oscillations are present; they
reduce at the last stages of the training. Also, not all parameters show the
same convergence rate since, for instance, the outlet number $3$ takes around
one-third of the total training epochs to converge to its reference value,
while the outlet number $5$ needs to approximately double the iterations to
have a same-accurate result.

Figure 6: Evolution of Parameter estimation statistics for each total resistances of the aorta in steady flow regime. All values were divided by their reference value and split vertically for better visualization. Dashed black lines represent the exact solution, and continuous color lines represent the estimation mean. The colored area surrounding the curves is determined by the minimum and maximum value found along the multiple realizations of the experiment |  Steady Problem |
---|---|---
|  $Q\ (\rm{mL/s})$ |  $P\ (\rm{mmHg})$ |  $R_{tot}\ (10^{3}\ \rm{dyn\cdot s\cdot cm^{-5}})$
Bnd. | Ref. | Mean | Std | Ref | Mean | Std | Ref | Mean | Std
$\Gamma_{in}$ | -298.88 | -297.16 | 2.07 | 248.39 | 248.92 | 0.63 | $-$ | $-$ | $-$
$\Gamma_{1}$ | 25.82 | 25.58 | 0.10 | 247.96 | 246.93 | 1.63 | 12.74 | 12.86 | 0.14
$\Gamma_{2}$ | 25.69 | 25.64 | 0.13 | 247.90 | 247.34 | 0.80 | 12.74 | 12.89 | 0.06
$\Gamma_{3}$ | 30.56 | 30.71 | 0.27 | 248.20 | 247.98 | 0.66 | 10.75 | 10.78 | 0.08
$\Gamma_{4}$ | 26.65 | 26.08 | 0.12 | 248.22 | 249.08 | 0.66 | 12.30 | 12.71 | 0.10
$\Gamma_{5}$ | 188.85 | 185.97 | 0.84 | 248.15 | 248.30 | 0.66 | 1.74 | 1.78 | 0.04
Table 4: Flow and mean pressure for each outlet of the aorta. The mean
pressures were computed as $P_{i}=\frac{1}{Area_{i}}\int_{\Gamma_{i}}pds$,
where $i$ is the corresponding outlet number. The mean and standard deviation
values reported are computed from the results of every seed.

Table 4 shows the final values of the estimation, the mean estimated flow, and
pressure at every outlet. These quantities were computed from the mean of the
velocities and pressures by every experiment realization. We have found that
the PINNs represent the solution at the inlet remarkably well, even though no
information was given at that boundary. Also, the estimated parameters match
their corresponding reference values well despite the lack of volumetric data
and the artifacts introduced to the slice measurements due to the added noise
and spatiotemporal interpolation. Finally, Figure 7 shows the velocity
streamlines computed from the inlet nodes of the mesh of both reference and
mean estimation velocities. This figure shows that the PINNs’ solution is
similar to the reference velocity while in the vessel’s lumen and moderately
fails for the nodes closer to the walls.

(a) Reference Velocity

(b) Mean Estimated Velocity

Figure 7: Velocity streamlines when the steady-state flow regime is imposed
for (a) the reference solution and (b) the mean estimated velocity represented
by the PINNs

###  3.2 Transient problem

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

Figure 8: Evolution of the parameter estimation statistics for each Windkessel outlet of the aorta in transient flow regime. All values were divided by their reference value and split vertically for better visualization. Dashed black lines represent the exact solution, and continuous color lines represent the estimation mean. The colored area surrounding the curves is determined by the minimum and maximum values found along the multiple realizations of the experiment. |  Transient Problem |
---|---|---
|  $R_{p}\ (10^{3}\ \rm{dyn\cdot s\cdot cm^{-5}})$ |  $R_{d}\ (10^{3}\ \rm{dyn\cdot s\cdot\rm{cm^{-5}}})$ |  $C_{d}\ (10^{-5}\ \rm{dyn^{-1}\cdot cm^{5}})$
Bnd. | Ref. | Mean | Std | Ref | Mean | Std | Ref | Mean | Std
$\Gamma_{1}$ | 0.71 | 1.02 | 0.10 | 12.02 | 15.76 | 1.87 | 8.26 | 6.42 | 0.70
$\Gamma_{2}$ | 0.71 | 1.10 | 0.12 | 12.02 | 16.75 | 1.58 | 8.26 | 6.01 | 0.52
$\Gamma_{3}$ | 0.60 | 0.91 | 0.15 | 10.14 | 14.51 | 1.88 | 9.78 | 6.99 | 0.82
$\Gamma_{4}$ | 0.69 | 0.84 | 0.14 | 11.61 | 14.46 | 1.94 | 8.55 | 7.01 | 0.80
$\Gamma_{5}$ | 0.10 | 0.12 | 0.02 | 1.65 | 2.02 | 0.22 | 60.15 | 50.18 | 5.51
Table 5: Parameter estimation obtained on the transient-flow regime.

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

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

tend to be overestimated in the upper aorta and underestimated in the downward
portion of the domain. Moreover, an oscillating behavior was found in the
flows and pressure in the outlets $\Gamma_{3}$ and $\Gamma_{4}$, outlets where
the standard deviation of the estimated parameters is the most. A FEM
simulation was also performed after the training, using the estimated
parameters and inlet velocity profile found in the experiment, which present
the median error among all realizations. This solution could be seen as a
post-process step of the PINNs’ estimation and generally presents reduced
oscillations compared with the raw estimated velocity and pressure fields. The
obtained curves of the post-process FEM simulation are depicted in colored
dashed lines.

Figure 9: Flows (red) and mean pressure (sky blue) curves at every aortic
outlet over time. Continuous lines represent mean estimation, while the area
plot is defined by the minimum and maximum value encountered in every
experiment realization. Black dashed lines correspond to the respective
reference value, while colored dashed lines correspond to the curves obtained
from the post-process FEM simulation.

Finally, streamlines computed from the nodes of $\Gamma_{in}$ were calculated
at two instants of the cardiac cycle and shown in Figure 10. This was done for
the reference velocity, the PINN’s estimation, and the FEM post-process
simulated velocity. Generally, the neural network can capture the main flow
features at systole but not so well during diastole, where the velocities are
lower and the measurements’ velocity-to-noise ratio is worse. When comparing
the reference velocity against the FEM post-process solution, no significant
gain in quality can be seen while in systole. However, during diastole, the
FEM solution captured the main flow features that the PINNs’ estimation failed
to represent, especially the vortices at the ascending aorta and the near-wall
velocity in the entire domain.

(a) Reference ($t=0.12$ s)

(b) PINNs ($t=0.12$ s)

(c) FEM with estimated
params. ($t=0.12$ s)

(d) Reference ($t=0.45$ s)

(e) PINNs ($t=0.45$ s)

(f) FEM with estimated
params. ($t=0.45$ s)

Figure 10: Velocity streamlines of the transient problem for the reference
solution (a and d), the mean estimated velocity estimated by the PINNs (b and
e), and the velocity obtained after a FEM simulation using the found
parameters (c and f) at two time instants: $t=0.12\ s$, which correspond to
peak systole, and at $t=0.45\ s$, during mid diastole.

##  4 Discussion

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

The stationary problem could be seen as a inexpensive version of the transient
problem, and in clinical hemodynamics, is usually studied when the transient
problem’s computational time exceeds the allowed limit imposed by the clinical
application itself [51, 52]. In terms of the quality of the estimations, it
could be seen from the results in Table 4 that although all parameters were
overestimated, the obtained values present a mean error of less than $2\%$
when comparing with their reference. On the other hand, when looking into the
reconstructed velocity streamlines in Figure 7, both reconstructed and
estimated fields are very similar, with some mild differences at the vessel
wall. One of the reasons for that difference is because FEM methods exactly
imposes the boundary condition on the walls, being the velocity there set
exactly to zero, as it happened with the reference solution. However, in the
case of the PINNs method presented here, the non-slip condition enters the
optimization as a loss term and, consequently, is not completely satisfied.
This produces some of the velocity lines to go out of the vessel through the
walls. This effect is also present in the transient flow regime.

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

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

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

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

##  5 Conclusions

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

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

## Acknowledgements

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

## References

  * [1] R. Erbel, F. Alfonso, C. Boileau, O. Dirsch, B. Eber, A. Haverich, H. Rakowski, J. Struyven, K. Radegran, U. Sechtem, et al., Diagnosis and management of aortic dissection: task force on aortic dissection, european society of cardiology, European Heart Journal 22 (18) (2001) 1642–1681.
  * [2] B. A. Carabello, W. J. Paulus, Aortic stenosis, The lancet 373 (9667) (2009) 956–966.
  * [3] N. Sakalihasan, R. Limet, O. D. Defawe, Abdominal aortic aneurysm, The Lancet 365 (9470) (2005) 1577–1589.
  * [4] D. W. McRobbie, MRI from Picture to Proton, Cambridge University Press, 2003.
  * [5] D. Fleischmann, A. S. Chin, L. Molvin, J. Wang, R. Hallett, Computed tomography angiography: a review and technical update, Radiologic Clinics 54 (1) (2016) 1–12.
  * [6] J. J. Hom, K. Ordovas, G. P. Reddy, Velocity-encoded cine MR imaging in aortic coarctation: functional assessment of hemodynamic events, Radiographics 28 (2) (2008) 407–416.
  * [7] D. Nolte, C. Bertoglio, Inverse problems in blood flow modeling: A review., International Journal for Numerical Methods in Biomedical Engineering 38 (8) (2022) e3613.
  * [8] A. Caiazzo, I. E. Vignon-Clementel, Mathematical modeling of blood flow in the cardiovascular system, Quantification of Biophysical Parameters in Medical Imaging (2018) 45–70.
  * [9] H. J. Kim, I. Vignon-Clementel, J. Coogan, C. Figueroa, K. Jansen, C. Taylor, Patient-specific modeling of blood flow and pressure in human coronary arteries, Annals of biomedical engineering 38 (10) (2010) 3195–3209.
  * [10] N. Grande Gutiérrez, T. Sinno, S. L. Diamond, A 1d–3d hybrid model of patient-specific coronary hemodynamics, Cardiovascular Engineering and Technology 1–12.
  * [11] C. Bertoglio, C. Conca, D. Nolte, G. Panasenko, K. Pileckas, Junction of models of different dimension for flows in tube structures by womersley-type interface conditions, SIAM Journal on Applied Mathematics 79 (3) (2019) 959–985.
  * [12] C. Bertoglio, A. Caiazzo, M. A. Fernández, Fractional-step schemes for the coupling of distributed and lumped models in hemodynamics, SIAM Journal on Scientific Computing 35 (3) (2013) B551–B575.
  * [13] M. Ismail, W. A. Wall, M. W. Gee, Adjoint-based inverse analysis of windkessel parameters for patient-specific vascular models, Journal of Computational Physics 244 (2013) 113–130.
  * [14] F. Marcinno’, A. Zingaro, I. Fumagalli, L. Dede’, C. Vergara, A computational study of blood flow dynamics in the pulmonary arteries, Vietnam Journal of Mathematics 51 (1) (2023) 127–149.
  * [15] R. Piersanti, F. Regazzoni, M. Salvador, A. F. Corno, C. Vergara, A. Quarteroni, et al., 3D–0D closed-loop model for the simulation of cardiac biventricular electromechanics, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114607.
  * [16] G. Pase, E. Brinkhuis, T. De Vries, J. Kosinka, T. Willems, C. Bertoglio, A parametric geometry model of the aortic valve for subject-specific blood flow simulations using a resistive approach, Biomechanics and Modeling in Mechanobiology 22 (3) (2023) 987–1002.
  * [17] N. Westerhof, J. Lankhaar, B. Westerhof, The arterial windkessel, Medical biological engineering and computing 47 (2009) 131–41.
  * [18] O. Franck, The basic shape of the arterial pulse. first treatise: mathematical analysis. 1899, J. Mol. Cell. Cardiol 22 (3) (1990) 255–277.
  * [19] C. J. Arthurs, N. Xiao, P. Moireau, T. Schaeffter, C. A. Figueroa, A flexible framework for sequential estimation of model parameters in computational hemodynamics, Advanced modeling and simulation in engineering sciences 7 (1) (2020) 1–37.
  * [20] E. Fevola, F. Ballarin, L. Jiménez-Juan, S. Fremes, S. Grivet-Talocia, G. Rozza, P. Triverio, An optimal control approach to determine resistance-type boundary conditions from in-vivo data for cardiovascular simulations, International Journal for Numerical Methods in Biomedical Engineering 37 (10) (2021) e3516.
  * [21] C. Bertoglio, P. Moireau, J.-F. Gerbeau, Sequential parameter estimation for fluid–structure problems: application to hemodynamics, International Journal for Numerical Methods in Biomedical Engineering 28 (4) (2012) 434–455.
  * [22] J. Garay, D. Nolte, M. Löcke, C. Bertoglio, Parameter estimation in fluid flow models from aliased velocity measurements, Inverse Problems 38 (9) (2022) 095002.
  * [23] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.
  * [24] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440.
  * [25] S. Cai, Z. Mao, Z. Wang, M. Yin, G. E. Karniadakis, Physics-informed neural networks (pinns) for fluid mechanics: A review., Acta Mechanica Sinica 37 (12) (2021) 1727–1738.
  * [26] H. Xiao, J.-L. Wu, S. Laizet, L. Duan, Flows over periodic hills of parameterized geometries: A dataset for data-driven turbulence modeling from direct simulations, Computers & Fluids 200 (2020) 104431.
  * [27] A. Mathews, M. Francisquez, J. W. Hughes, D. R. Hatch, B. Zhu, B. N. Rogers, Uncovering turbulent plasma dynamics via deep learning from partial observations, Physical Review E 104 (2) (2021) 025205.
  * [28] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, E. Kuhl, Physics-informed neural networks for cardiac activation mapping, Frontiers in Physics 8 (2020) 42.
  * [29] C. Ruiz Herrera, T. Grandits, G. Plank, P. Perdikaris, F. Sahli Costabal, S. Pezzuto, Physics-informed neural networks to learn cardiac fiber orientation from multiple electroanatomical maps, Engineering with Computers 38 (5) (2022) 3957–3973.
  * [30] W. Li, M. Z. Bazant, J. Zhu, A physics-guided neural network framework for elastic plates: Comparison of governing equations-based and energy-based approaches, Computer Methods in Applied Mechanics and Engineering 383 (2021) 113933\.
  * [31] S. Maharjan, B. Guidio, A. Fathi, C. Jeong, Deep and convolutional neural networks for identifying vertically-propagating incoming seismic wave motion into a heterogeneous, damped soil column, Soil Dynamics and Earthquake Engineering 162 (2022) 107510.
  * [32] V. Tac, K. Linka, F. Sahli Costabal, E. Kuhl, A. B. Tepole, Benchmarks for physics-informed data-driven hyperelasticity, arXiv preprint arXiv:2301.10714.
  * [33] M. Bazmara, M. Silani, M. Mianroodi, M. Sheibanian, Physics-informed neural networks for nonlinear bending of 3D functionally graded beam, Structures 49 (2023) 152–162.
  * [34] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A deep learning framework for solution and discovery in solid mechanics, arXiv preprint arXiv:2003.02751.
  * [35] A. Arzani, J.-X. Wang, R. M. D’Souza, Uncovering near-wall blood flow from sparse data with physics-informed neural networks, Physics of Fluids 33 (7) (2021) 071905.
  * [36] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4D flow MRI data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112623\.
  * [37] A. Arzani, J.-X. Wang, M. S. Sacks, S. C. Shadden, Machine learning for cardiovascular biomechanics modeling: Challenges and beyond, Annals of Biomedical Engineering (2022) 1–13.
  * [38] L. Sun, H. Gao, S. Pan, J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Computer Methods in Applied Mechanics and Engineering 361 (2020) 112732.
  * [39] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through physics–informed neural networks: where we are and what’s next, Journal of Scientific Computing 92 (3) (2022) 88.
  * [40] F. Saul, D. Klaus, Y. Aristidou, A. Wiemeyer, B. Lösse, Non-invasive oscillometric wrist and upper arm blood pressure measurements compared with invasive values., Zeitschrift fur Kardiologie 85 (1996) 127–129.
  * [41] H.-M. Cheng, K.-L. Wang, Y.-H. Chen, S.-J. Lin, L.-C. Chen, S.-H. Sung, P. Y.-A. Ding, W.-C. Yu, J.-W. Chen, C.-H. Chen, Estimation of central systolic blood pressure using an oscillometric blood pressure monitor, Hypertension Research 33 (6) (2010) 592–599.
  * [42] M. Forouzanfar, H. R. Dajani, V. Z. Groza, M. Bolic, S. Rajan, I. Batkin, Oscillometric blood pressure estimation: past, present, and future, IEEE reviews in biomedical engineering 8 (2015) 44–63.
  * [43] K. Jänich, L. Kay, Vector analysis, Vol. 5, Springer, 2001.
  * [44] S. Wang, Y. Teng, P. Perdikaris, Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal on Scientific Computing 43 (5) (2021) A3055–A3081.
  * [45] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980.
  * [46] P. Ramachandran, B. Zoph, Q. V. Le, Searching for activation functions, arXiv preprint arXiv:1710.05941.
  * [47] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in: Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.
  * [48] N. M. Wilson, A. K. Ortiz, A. B. Johnson, The vascular model repository: a public resource of medical imaging data and blood flow simulation results, Journal of Medical Devices 7 (4).
  * [49] Y. Bazilevs, J. Gohean, T. Hughes, R. Moser, Y. Zhang, Patient-specific isogeometric fluid–structure interaction analysis of thoracic aortic blood flow due to implantation of the jarvik 2000 left ventricular assist device, Computer Methods in Applied Mechanics and Engineering 198 (45-46) (2009) 3534–3550.
  * [50] P. Irarrazaval, A. Dehghan Firoozabadi, S. Uribe, C. Tejos, C. Sing-Long, Noise estimation for the velocity in mri phase-contrast, Magnetic Resonance Imaging 63 (2019) 250–257.
  * [51] R. Vinoth, D. Kumar, R. Adhikari, S. Vijayapradeep, K. Geetha, R. Ilavarasi, S. Mahalingam, Steady and transient flow cfd simulations in an aorta model of normal and aortic aneurysm subjects, in: The Proceedings of the International Conference on Sensing and Imaging, Springer, 2019, pp. 29–43.
  * [52] F. Ene, P. Delassus, L. Morris, The influence of computational assumptions on analysing abdominal aortic aneurysm haemodynamics, Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 228 (8) (2014) 768–780.
  * [53] H. Gao, L. Sun, J.-X. Wang, Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels, Physics of Fluids 33 (7).
  * [54] D. Nolte, C. Bertoglio, Reducing the impact of geometric errors in flow computations using velocity measurements, International journal for numerical methods in biomedical engineering 35 (6) (2019) e3203.
  * [55] A. T. Mohan, N. Lubbers, D. Livescu, M. Chertkov, Embedding hard physical constraints in neural network coarse-graining of 3d turbulence, arXiv preprint arXiv:2002.00021.
  * [56] F. Marner, et al., Potential-based formulations of the navier-stokes equations and their application, Ph.D. thesis, Durham University (2019).
  * [57] M. F. Fathi, I. Perez-Raya, A. Baghaie, P. Berg, G. Janiga, A. Arzani, R. M. D’Souza, Super-resolution and denoising of 4D-flow MRI using physics-informed deep neural nets, Computer Methods and Programs in Biomedicine 197 (2020) 105729.
  * [58] J. Mura, A. M. Pino, J. Sotelo, I. Valverde, C. Tejos, M. E. Andia, P. Irarrazaval, S. Uribe, Enhancing the velocity data from 4d flow mr images by reducing its divergence, IEEE transactions on medical imaging 35 (10) (2016) 2353–2364.
  * [59] J. Garay, H. Mella, J. Sotelo, C. Cárcamo, S. Uribe, C. Bertoglio, J. Mura, A new mathematical model for verifying the navier-stokes compatibility of 4d flow mri data.
  * [60] L. Partin, D. E. Schiavazzi, C. A. S. Long, An analysis of reconstruction noise from undersampled 4d flow mri, Biomedical Signal Processing and Control 84 (2023) 104800.