# BP-DeepONet: A new method for cuffless blood pressure estimation using the
physcis-informed DeepONet

Lingfeng Li<EMAIL_ADDRESS>Xue-Cheng Tai<EMAIL_ADDRESS>Raymond
Hon-Fu Chan<EMAIL_ADDRESS>
###### Abstract

Cardiovascular diseases (CVDs) are the leading cause of death worldwide, with
blood pressure serving as a crucial indicator. Arterial blood pressure (ABP)
waveforms provide continuous pressure measurements throughout the cardiac
cycle and offer valuable diagnostic insights. Consequently, there is a
significant demand for non-invasive and cuff-less methods to measure ABP
waveforms continuously. Accurate prediction of ABP waveforms can also improve
the estimation of mean blood pressure, an essential cardiovascular health
characteristic. This study proposes a novel framework based on the physics-
informed DeepONet approach to predict ABP waveforms. Unlike previous methods,
our approach requires the predicted ABP waveforms to satisfy the Navier-Stokes
equation with a time-periodic condition and a Windkessel boundary condition.
Notably, our framework is the first to predict ABP waveforms continuously,
both with location and time, within the part of the artery that is being
simulated. Furthermore, our method only requires ground truth data at the
outlet boundary and can handle periodic conditions with varying periods.
Incorporating the Windkessel boundary condition in our solution allows for
generating natural physical reflection waves, which closely resemble
measurements observed in real-world cases. Moreover, accurately estimating the
hyper-parameters in the Navier-Stokes equation for our simulations poses a
significant challenge. To overcome this obstacle, we introduce the concept of
meta-learning, enabling the neural networks to learn these parameters during
the training process. To demonstrate the effectiveness of our approach, we
conduct numerical experiments that showcase its superiority over traditional
methods, which typically only predict systolic blood pressure (SBP) and
diastolic blood pressure (DBP). Our proposed method not only accurately
predicts continuous ABP waveforms but also outperforms these traditional
approaches.

††journal: xxxxx

[1]organization=Hong Kong Centre for Cerebro-cardiovascular Health
Engineering, city=Hong Kong, country=China [2]organization=Norwegian Research
Center, city=Bergen, country=Norway [3]organization=Department of mathematics,
City university of Hong Kong, city=Hong Kong, country=China

##  1 Introduction

According to the World Health Organization (WHO), cardiovascular diseases
(CVDs) are the leading cause of death globally [1]. In 2019, there were about
32% of all global deaths caused by CVDs. The majority of these deaths were due
to heart attack and stroke. There is sufficient evidence to support that high
blood pressure (or hypertension) is the predominant cause of CVDs [2]. In
practice, the morphology of the arterial blood pressure (ABP) waveforms can
reflect the cardiovascular status of humans [3, 4, 5, 6]. Therefore, it is of
great importance to monitor ABP waveforms continuously to detect and prevent
diseases.

In the clinical setting, the gold standard for ABP waveform monitoring is the
invasive method, which involves implanting an invasive pressure sensor into
the artery. The invasive method is very risky, and the used devices are
expensive. It is commonly used in intensive care units (ICUs). On the other
hand, the non-invasive method is a more desirable choice for the daily
monitoring of ABP waveforms.

The deep learning-based method for non-invasive ABP waveform estimation has
attracted huge research interest in recent years. These methods usually
estimate the ABP from other physiological signals that are much easier to
measure. Electrocardiogram (ECG) and photoplethysmogram (PPG) are the most
commonly used signals. The ECG signal records the heart’s electrical activity
through cardiac cycles. It is measured by some electrodes, which can detect
the small electrical changes resulting from the cardiac activity. The patterns
of ECG signals can indicate many cardiac abnormalities. The PPG signal is an
optical signal that can detect blood volume changes in the microvascular bed
of tissue. It is often obtained by illuminating the skin with a pulse oximeter
and measuring changes in light absorption. Examples of filtered ECG and PPG
signals are shown in Figure 1.

Figure 1: Filtered ECG and PPG signals in 10 seconds.

Deep learning models were introduced to estimate systolic and diastolic blood
pressure in [7]. Then, more complicated neural networks are used to estimate
the entire ABP waveform [8, 9, 10]. In [8], the authors proposed a neural
network consisting of a convolutional neural network (CNN), a bi-directional
gated recurrent unit (Bi-GRU), and attention layers to estimate ABP from ECG
and PPG. In [9], a network that predicts ABP only using PPG signals was
proposed. In [10], the authors introduced a modified transformer-based model
to predict the ABP waveform from ECG, PPG, and some morphological extracted
features.

The models above are all trained in a data-driven way, i.e., trained by
minimizing the difference between the ground truth label and the network
prediction. Since the labeled data is usually measured by invasive methods in
the radial artery [11], the trained neural network will only predict the
radial ABP waveform. However, ABP waveforms differ at different locations in
the arterial systems [12, 13] (See Figure 2), and radial waveforms may not be
sufficient to diagnose CVD and reflect the status of the cardiovascular
system. For example, some studies have shown that the brachial waveform is
more valuable than the radial waveform when assessing cardiac contractility
[14, 15]. Another important observation is the pulse pressure amplification
effect. Pulse pressure is defined as the difference between the maximum and
minimum of the pressure waveforms, and this difference will increase through
the propagation of the pulse wave in the arterial system. All current deep
learning-based algorithms can only predict radial ABP waveforms since they are
trained with labeled data measured at the radial artery.

Figure 2: Illustration of the pulse pressure amplification in arterial
systems.

In this work, we aim to develop new methods to predict ABP waveforms at
different locations and at all times in the artery conditioned to the physical
law that governs the blood flow, i.e. the Navier-Stokes equation. We use the
physics-informed DeepONet (PI-DeepONet) [16]. DeepONet was initially proposed
to learn the solution operator that maps variable inputs to the corresponding
latent solutions of certain partial differential equation (PDE) systems [17],
and its training process is entirely data-driven. If the underlying PDE system
is explicitly given, then the PDE can be incorporated into the training loss
like the physics-informed neural networks (PINN) [18]. The DeepONet trained in
this way is called the physics-informed DeepONet. In the arterial system, the
hemodynamics, i.e., blood flow dynamics, can be modeled and simulated by the
Navier-Stokes equation [19, 20, 21]. Therefore, we can incorporate the Navier-
Stokes equation into the training loss and learn the solution operator that
maps physiological signals to solutions of the corresponding PDEs. In our
setting, the spatial domain of the Navier-Stokes equation is a segment of the
artery, and the labeled data can serve as a Dirichlet boundary condition.
After providing other necessary initial/boundary conditions, we can identify
the PDE solution in the whole domain. Consequently, we can obtain the ABP
waveforms at different locations and at all times.

We summarize our main contributions as follows:

  * 1.

This work is the first attempt to apply physics-informed operator learning to
ABP waveform estimation. The proposed methods are the first ones to predict
ABP waveforms at different locations in the artery conditioned by the physical
law that governs the blood flow, i.e., the Navier-Stokes equation. The
proposed methods enforce neural network solutions to satisfy a 1-D Navier-
Stokes equation with a Windkessel boundary condition and a time-periodic
condition. The Windkessel boundary condition can simulate physical reflection
waves.

  * 2.

The proposed methods can predict the blood pressure waveforms with reasonably
good accuracy and preserve some physical properties like pulse pressure
amplification.

  * 3.

To the best of our knowledge, the proposed methods are the first ones that
incorporate variable time-periodic conditions to DeepONet. We also propose an
effective implementation method for the network training.

  * 4.

We also incorporate meta-learning into the training to estimate model hyper-
parameters automatically from the input physiological signals.

  * 5.

During training, we only require the blood pressure measurement at the outlet
side of the domain, and no extra measurement is needed.

##  2 One-dimensional model for blood flow simulation

The hemodynamics in a vessel segment can be effectively modeled and simulated
by the one-dimensional Navier-Stokes equation [20, 19]. The one-dimensional
model is a simplification of the three-dimensional Navier-Stokes equation by
making several assumptions on the computational domain:

  * 1.

The vessel is a straight cylinder with the axis oriented along the coordinate
$z$ direction (See Figure 3).

  * 2.

The vessel wall displaces along the radial direction only, which implies that
each axial cross-section is always circular.

  * 3.

The blood pressure and blood flow rate are constant within each axial cross-
section.

  * 4.

All body forces, like the gravity, are neglected.

Figure 3: Simplification of the three-dimensional model to one-dimensional
model. The vessel is assumed to be a straight cylinder, and each axial cross-
section is circular.

Thus, the hemodynamics can be described by three quantities, namely the blood
flow rate $Q(z,t)$, the axial cross-sectional area $A(z,t)$, and the blood
pressure $P(z,t)$, where $z\in[0,L]$ is the spatial coordinates and
$t\in[0,T]$ is the temporal coordinates. $L$ denotes the length of the vessel
segment, and $T$ represents the time duration. Then, the one-dimensional model
can be written as:

$\displaystyle\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,$
(1) $\displaystyle\frac{\partial Q}{\partial t}+\frac{\partial}{\partial
z}(\frac{Q^{2}}{A})+\frac{A}{\rho}\frac{\partial P}{\partial
z}+K_{r}\frac{Q}{A}=0,$ (2)

where $\rho$ is the blood density, and $K_{r}$ is a resistance parameter
related to blood viscosity. To close this system, the blood pressure $P$ and
the cross-sectional area are assumed to satisfy:

$P=P_{ext}+\beta\frac{\sqrt{A}-\sqrt{A_{0}}}{A_{0}},$ (3)

where $P_{ext}$ is the external pressure, $A_{0}$ is the cross-sectional area
at the reference state, and $\beta$ is a constant depending on the physical
and mechanical properties of vessels. More details of the model derivation can
be found in [20, 19]. A simple characteristic analysis shows that the forward
characteristic variable $W_{1}$ and backward characteristic variable $W_{2}$
in the system (1)-(2) are:

$W_{1}(A,Q)=\frac{Q}{A}+4\sqrt{\frac{\beta}{2\rho
A_{0}}}(A^{1/4}-A_{0}^{1/4}),$

and

$W_{2}(A,Q)=\frac{Q}{A}-4\sqrt{\frac{\beta}{2\rho
A_{0}}}(A^{1/4}-A_{0}^{1/4}).$

The forward wave and backward wave are traveling at the speed

$\lambda_{1}(A,Q)=\frac{Q}{A}+\sqrt{\frac{\beta}{2\rho A_{0}}}A^{1/4}$

and

$\lambda_{2}(A,Q)=\frac{Q}{A}-\sqrt{\frac{\beta}{2\rho A_{0}}}A^{1/4}$

respectively. In practice, we have

$\sqrt{\frac{\beta}{2\rho A_{0}}}A^{1/4}\gg Q/A,$

so $W_{1}$ and $W_{2}$ always travel in opposite directions. The speed of the
forward wave $\lambda_{1}\approx\sqrt{\frac{\beta}{2\rho A_{0}}}A^{1/4}$ is
often called pulse wave velocity.

Using the relation (3), we can rewrite the equations (1)-(2) as

$\frac{\partial U}{\partial t}+H(U)\frac{\partial U}{\partial
z}=B(U),\quad(z,t)\in[0,L]\times[0,T]$ (4)

where

$U=\begin{pmatrix}P\\\ Q\end{pmatrix},$
$H(U)=\begin{pmatrix}0&\frac{\beta}{2a(P)A_{0}}\\\
\frac{a(P)^{2}}{\rho}-\frac{2A_{0}Q^{2}}{\beta
a^{3}(P)}&\frac{2Q}{a(P)^{2}}\end{pmatrix},$ $B(U)=\begin{pmatrix}0\\\
-K_{r}\frac{Q}{a(P)^{2}}\end{pmatrix},$

and

$a(P)=\frac{A_{0}}{\beta}(P-P_{ext})+\sqrt{A_{0}}.$

Even though the one-dimensional model can not capture the blood flow details,
it can effectively describe the wave propagation within the vessel. This model
has been validated by both in vitro data [22, 23] and in vivo data [24, 25].
There are many effective numerical solvers, like the MacCormack scheme[26,
27], Taylor-Galerkin scheme [28, 29], MUSCL [30, 31], and the local
discontinuous Galerkin scheme [32, 33, 34]. A comparison of different
numerical methods has also been done in [35].

##  3 Physics-informed machine learning for PDEs

Machine learning algorithms for PDEs have been extensively studied in recent
years. In general, there are two main directions in this research field:
learning an instance of PDEs and learning the solution operator of PDEs. The
most famous work along the first direction is the physics-informed neural
networks (PINNs) proposed in [18]. This work aims to solve an instance of PDEs
using a neural network. Compared to traditional numerical solvers, the PINN is
mesh-free, easy to implement, and can solve high-dimensional problems.
However, it can not achieve very high accuracy for non-linear PDEs. The second
direction aims to learn the solution operator that maps from some input spaces
to the solution space using some neural networks. This type of network is also
called the neural operator. The most popular neural operators include the
DeepONet [17] and the Fourier Neural operator (FNO) [36]. The learned network
can quickly predict the PDE solution corresponding to the given input. We will
provide a general introduction to two types of methods in this section.

###  3.1 Physics-informed neural networks (PINNs)

Let’s consider a general form of a $d$-dimensional boundary value problem:

$\displaystyle\mathcal{L}(u)(x)$ $\displaystyle=f(x),\quad x\in\Omega$
$\displaystyle\mathcal{B}(u)(x)$ $\displaystyle=g(x),\quad x\in\partial\Omega$

where $\Omega\in\mathbb{R}^{d}$ is the problem domain, $\mathcal{L}(u)$
denotes a general differential operator, and $\mathcal{B}(u)$ denotes
corresponding boundary conditions. The PINN method first defines a neural
network to approximate the solution to this specific instance of PDE:
$\mathcal{N}(x;\theta):\Omega\rightarrow\mathbb{R}$, where $\theta$ denotes
the set of all trainable parameters in the network. Then, this network is
trained by solving the optimization problem:

$\min_{\theta}\sum_{j=1}^{N_{0}}|\mathcal{L}(\mathcal{N}(x_{j};\theta))-f(x_{j})|^{2}+\omega\sum_{j=1}^{N_{1}}|\mathcal{B}(\hat{x}_{j})-g(\hat{x}_{j})|^{2}$

where $\omega$ is the weight of the boundary residual,
$\\{x_{j}\\}_{j=1}^{N_{0}}$ is collocation points randomly sampled from
$\Omega$, and $\\{\hat{x}_{j}\\}_{j=1}^{N_{1}}$ is collocation points randomly
generated from $\partial\Omega$. The partial derivatives of $\mathcal{N}$ can
be easily calculated through the back-propagation mechanism. An efficient
implementation is provided in the DeepXDE package [37] for both PyTorch and
TensorFlow platforms. Some similar algorithms in the literature use similar
ideas, such as the deep Ritz method (DRM) [38] and the mixed residual method
(MIM) [39, 40]. The DRM uses the variational formulation of PDEs to train
networks, and the MIM transforms high-order PDEs into first-order systems to
train networks. Various numerical analyses have also been conducted for these
physics-informed methods [41, 42, 43, 44, 45].

One popular network structure for the PINN method is the fully connected
network (FCN) (Figure 4(a)). An FCN is typically defined as compositions of
many fully connected layers:

$\displaystyle\mathcal{N}_{FCN}(x;\theta):=(\mathcal{F}_{L}\circ\mathcal{F}_{L-1}\circ\dots\circ\mathcal{F}_{1})(x),$
$\displaystyle\mathcal{F}_{l}(z):=\sigma_{l}(W_{l}z+b_{l}),\quad\forall
l=1,\dots,L-1,$ $\displaystyle W_{l}\in\mathbb{R}^{d_{l,in}\times
d_{l,out}},\quad b_{l}\in\mathbb{R}^{d_{l,out}},$ $\displaystyle
d_{1,in}=d,\quad d_{l,in}=d_{l-1,out},\quad l=2,\dots,L-1,$
$\displaystyle\mathcal{F}_{L}(z):=\sigma_{L}(W_{L}z+b_{L}),\quad
W_{L}\in\mathbb{R}^{d_{L-1,out},1},\quad b_{L}\in\mathbb{R},$

where $\sigma_{l}$ is the activation function for each layer, $d_{l,in}$
denotes the input dimension of the $l$th layer, and $d_{l,out}$ denotes the
output dimension of the $l$th layer. Choices for the activation function
usually include Sigmoid, Tanh, ReLU, etc. For the last layer $F_{L}$, the
activation function can be the identify function, i.e., $\sigma_{L}(z)=z$.
When increasing the depth $L$, the training of the FCN may become slower
because of the vanishing gradient effect. To overcome this limitation, we can
use the residual network incorporating the skip connection structure [46]. A
residual network (Figure 4(b)) is defined as:

$\displaystyle\mathcal{N}_{Res}(x;\theta):=(\hat{\mathcal{F}}_{L}\circ\hat{\mathcal{F}}_{L-1}\circ\dots\circ\hat{\mathcal{F}}_{1})(x),$
$\displaystyle\hat{\mathcal{F}}_{l}(z):=\sigma_{l}(W_{l,2}\sigma_{l}(W_{l,1}z+b_{l,1})+b_{l,2}+z),\quad\forall
l=1,\dots,L-1,$ $\displaystyle W_{l,1}\in\mathbb{R}^{d_{l,in}\times
d_{l,out}},\quad W_{l,2}\in\mathbb{R}^{d_{l,out}\times d_{l,out}},\quad
b_{l,1},b_{l,2}\in\mathbb{R}^{d_{l,out}},$ $\displaystyle d_{1,in}=d,\quad
d_{l,in}=d_{l-1,out},\quad l=2,\dots,L-1,$
$\displaystyle\mathcal{F}_{L}(z):=\sigma_{L}(W_{L}z+b_{L}),\quad
W_{L}\in\mathbb{R}^{d_{L-1,out},1},\quad b_{L}\in\mathbb{R}.$

(a) Fully connected network

(b) Residual network

Figure 4: Structures of the fully connected network and the residual network.

###  3.2 Neural operators for PDEs

Let’s consider a class of parametric PDEs:

$\displaystyle\mathcal{L}(u)(x;\alpha)$ $\displaystyle=f(x;\alpha),\quad
x\in\Omega$ $\displaystyle\mathcal{B}(u;\alpha)(x)$
$\displaystyle=g(x;\alpha),\quad x\in\partial\Omega,$

where $\alpha$ is the parameters that define the problem; for example,
$\alpha$ can be the diffusion coefficient in a convection-diffusion equation.
The neural operator aims to approximate the solution operator
$\mathcal{M}:\mathcal{A}\rightarrow\mathcal{U}$, where $\mathcal{A}$ is the
function space for $\alpha$, and $\mathcal{U}$ is the space for the solution
$u$. In practice, the input space can also be other function spaces, e.g., the
space of the source term $f(x)$. Two types of neural operators have been
extensively applied to various tasks in the literature: the DeepONet [17] and
the FNO [36].

A plain DeepONet (Figure 5(a)) consists of two sub-networks: a branch net and
a trunk net. It is defined as

$\mathcal{N}_{DON}(\hat{\alpha},x):=Branch(\hat{\alpha})^{\top}Trunk(x)$

where $\hat{\alpha}$ is the discretization of $\alpha$ on a given mesh, $x$ is
any point in $\Omega$. The outputs of the branch net and the trunk net are
$p$-dimensional vectors, where $p$ is a user-specified hyper-parameter. The
DeepONet $\mathcal{N}(\hat{\alpha},x)$ gives a prediction to
$\mathcal{M}(\alpha)(x)$. To predict the solution on multiple points, we need
to evaluate the trunk net multiple times.

A plain FNO (Figure 5(b)) is the composition of many Fourier layers:

$\displaystyle\mathcal{N}_{FNO}(\hat{\alpha}):=(\tilde{\mathcal{F}}_{L}\circ\tilde{\mathcal{F}}_{L-1}\dots\circ\tilde{\mathcal{F}}_{1})(\hat{\alpha})$
$\displaystyle\tilde{\mathcal{F}}_{l}(z):=\sigma_{l}(W_{l}z+FFT^{-1}(R_{l}(FFT(z))))$

where $\sigma_{l}$ is the activation function, $FFT$ is the fast Fourier
transformation, $FFT^{-1}$ is the inverse fast Fourier transformation, $W_{l}$
is a linear transformation, and $R_{l}(z)$ is an operator that includes a
truncation of first $m$ frequency modes and mode-wise complex matrix
multiplications. The output of the FNO is of the same size as $\hat{\alpha}$.

(a) DeepONet

(b) FNO

Figure 5: Structures of the DeepONet and FNO.

The training of neural operators is fully data-driven. Suppose we have a set
of $\\{\hat{\alpha}_{i}\\}_{i=1}^{N}$ from $\mathcal{A}$ and corresponding PDE
solutions $\\{\hat{u}_{i}\\}_{i=1}^{N}$ defined on a discretization
$\\{x_{j}\\}_{j=1}^{M}$. Then we train the neural operator by minimizing the
$L_{2}$ difference between the network outputs and the reference solution:

$\min_{\theta}\frac{1}{MN}\sum_{i=1}^{N}\sum_{j=1}^{M}|\mathcal{N}_{DON}(\hat{\alpha}_{i},x_{j})-\hat{u}_{i}(x_{j})|^{2}$

for the DeepONet, or

$\min_{\theta}\frac{1}{N}\sum_{i=1}^{N}|\mathcal{N}_{FNO}(\hat{\alpha}_{i})-\hat{u}_{i}|_{2}^{2}$

for the FNO. A comprehensive study of two different neural operators is given
in [47], which includes numerical experiments on various linear and non-linear
PDEs. In general, the performances of two types of neural operators are
comparable in various problems. The FNO is more efficient in training and
testing, while the DeepONet has a more flexible structure that can be easily
applied to different tasks.

Besides the data-driven training method, we can apply the physics-informed
loss like PINN to train neural operators. Such methods are called physics-
informed neural operators [16, 48, 49]. In many engineering problems,
obtaining the entire solution to PDEs may be very difficult or expensive.
Consequently, it is more desirable to incorporate the physics laws into the
training process to regularize the neural operators. We can construct a
similar residual loss as the PINN method for each pair of training samples and
minimize the summation over all samples. For the DeepONet, we can use back-
propagation to evaluate the derivatives. For the FNO, we can use numerical
differentiation to evaluate the derivatives.

##  4 The physics informed DeepONet for ABP waveform estimation

###  4.1 Problem setup

This work aims to learn the operator that maps the physiological signals to
the solution of the Navier-Stokes equation:

$\mathcal{O}:\mathcal{S}([0,T])\rightarrow\mathcal{U}([0,L]\times[0,T])$

where $\mathcal{S}([0,T])$ is the set of physiological signals of length $T$
and $\mathcal{U}([0,L]\times[0,T])$ is the space of solutions to Navier-Stokes
equations defined on $[0,L]\times[0,T]$.

Suppose we have a training dataset consisting of $N$ pairs of samples:

$\\{(s_{i},p_{i})\\}_{i=1}^{N},$

where $s_{i}$ is a set of physiological signals, and $p_{i}$ is the
corresponding ABP waveform measured in the radial artery that lies in the
forearm. All waveform data have the same length. The physiological signal
$s_{i}$ serves as the input to the neural network and the blood pressure
waveform $p_{i}$ is used to define the boundary condition.

For each sample $i=1,\dots,N$, its hemodynamics is defined by
$U_{i}=(P_{i},Q_{i})$. We assume each $U_{i}$ satisfies the Navier-stokes
system described in (4):

$\frac{\partial U_{i}}{\partial t}+H(U_{i})\frac{\partial U_{i}}{\partial
z}=B(U_{i}),\quad(z,t)\in[0,L]\times[0,T].$ (5)

We also need to provide suitable initial and boundary conditions to identify
unique PDE solutions. When solving this system using some traditional
numerical methods, we usually need to impose an inlet boundary condition and
an outlet boundary condition. The inlet boundary condition is usually set as a
Dirichlet boundary condition which prescribes the value of the blood pressure
$P$ or blood flow rate $Q$ at the inlet side. The outlet boundary condition
can determine the reflection waves that run backward in the artery. A simple
choice for the outlet boundary condition is the non-reflecting boundary
condition [20]: $W_{2}=0$ where $W_{2}$ is the backward characteristic
variable. Though the non-reflecting boundary condition is very easy to
implement, it can not provide a physically correct reflecting wave. A better
choice is the three-element Windkessel boundary condition [50]:

$W(P_{i},Q_{i}):=Q_{i}(L,t)(1+\frac{R_{1}}{R_{2}})+CR_{1}\frac{\partial
Q_{i}}{\partial t}(L,t)-\frac{P_{i}(L,t)}{R_{2}}-C\frac{\partial
P_{i}}{\partial t}(L,t)=0,\quad t\in[0,T].$ (6)

where $R_{1}$, $R_{2}$, and $C$ are some hyper-parameters related to the
artery properties.

######  Remark.

The Windkessel model describes the hemodynamics of the arterial system in
terms of resistance and compliance, and it is widely used as an outflow
boundary condition in the numerical simulation of human arterial systems [22,
51, 21]. The Navier-Stokes equation with this Windkessel boundary condition
can produce physical reflection waves in the solutions. A simple comparison of
two outflow boundary conditions is shown in Figure 6 where we can observe that
the blood pressure simulated using the non-reflecting boundary condition is
not physically correct.

(a) Prescribed inlet boundary condition

(b) Simulated blood pressure using the non-reflecting boundary condition

(c) Simulated blood pressure using the Windkessel boundary condition

Figure 6: Comparison of the simulated blood pressure using different outflow
boundary conditions. (a) is the inlet boundary condition given at $z=0$; (b)
is the simulated blood pressure at $z=L$ using the non-reflecting outflow
boundary condition; and (c) is the simulated blood pressure at $z=L$ using the
3-element Windkessel outflow boundary condition. (b) and (c) are computed
using a standard MacCormack scheme.

In our case, since the pressure data $p_{i}$, $i=1,\dots,N$, is measured near
the end of radial arteries, it is better to formulate it also as an outlet
boundary condition:

$P_{i}(L,t)=p_{i}(t),\quad t\in[0,T].$ (7)

To better simulate the reflection waves in the artery, we use the three-
element Windkessel model (6) as another boundary condition.

For the initial condition, we impose a time-periodic condition on both $P_{i}$
and $Q_{i}$ instead of giving a prescribed initial value:

$U_{i}(z,t+\delta_{i})=U_{i}(z,t),\quad(z,t)\in[0,L]\times[0,T]$ (8)

where $\delta_{i}$ is the period for the $i$-th training sample, and it can be
estimated from the physiological signals $s_{i}$. Typically, the waveforms are
almost periodic within a short time.

######  Remark.

Conventional methods directly prescribe the initial value of solutions. One
commonly used initial condition is $U(z,0)=(P_{ext},0)$. However, it is very
difficult to measure the value of $P_{i}$ and $Q_{i}$ at the initial time for
each sample and it is not realistic to assume an initial condition like
$U(z,0)=(P_{ext},0)$. Therefore, using a time-periodic condition is more
suitable in our situation.

By solving the system (4) with the boundary condition (7) (6) and initial
condition (8), we would be able to predict the hemodynamics in the radial
artery. For conventional numerical methods, it would be difficult to solve
under this setting, because there is no prescribed initial value given and
there is no boundary condition at the inlet side; but it is easy to handle by
the physics-informed method.

In this problem, the period for each sample is different, which will cause
some difficulty in computing the loss function. In this work, we will propose
a new efficient implementation method to handle this issue.

###  4.2 BP-DeepONet

We propose a new physics-informed DeepONet, BP-DeepONet, to learn the operator
(4.1). We adopt the same structure as the original DeepONet. For the branch
net, we use a structure that combines a one-dimensional convolutional ResNet
and a bi-directional LSTM, see Figure 7. The inputs to the branch net are some
segments of physiological signals $s_{i}$ and the output would be two latent
vectors $\mathbf{b_{1}}$ and $\mathbf{b_{2}}$. conv1d(channels=$c$,kernel=$s$)
represents a standard one-dimensional convolutional layer whose number of
output channels is $c$ and kernel size is $s$. Bi-LSTM (size=$h$,layers=$l$)
represents a bi-directional LSTM block whose input size and hidden size are
$h$, and the number of layers is $l$. A similar structure has already been
used in [52] to extract features from physiological signals to predict ABP
waveforms. For the trunk net, we adopt a simple, fully connected ResNet as
Figure 4(b) with the feature expansion (Figure 8). Feature expansion [53] is a
technique to impose the periodic condition on PINN. To impose the time-
periodic constraints (8), we can simply extend the inputs $(z,t)$ to

$(z,\sin(\frac{2\pi t}{\delta_{i}}),\cos(\frac{2\pi
t}{\delta_{i}}),\sin(\frac{4\pi t}{\delta_{i}}),\cos(\frac{4\pi
t}{\delta_{i}}),\dots).$

Then, our network outputs will be strictly time-periodic with period
$\delta_{i}$. The output of the trunk net is a vector $\mathbf{k}$ whose
dimension is the same as $\mathbf{b_{1}}$ and $\mathbf{b_{2}}$. The final
prediction of the BP-DeepONet will be calculated as follows:

$P_{\theta}(s_{i},\delta_{i},(z,t))=\mathbf{b}_{1}(s_{i})^{\top}\mathbf{k}(z,t,\delta_{i}),\quad
Q_{\theta}(s_{i},\delta_{i},(z,t))=\mathbf{b}_{2}(s_{i})^{\top}\mathbf{k}(z,t,\delta_{i})$

where $\theta$ is the set of all trainable parameters in the BP-DeepONet.

Figure 7: Structure of the branch net. Figure 8: Structure of the trunk net.

The training loss function of our network is defined as

$\mathcal{L}(P_{\theta},Q_{\theta})=\mathcal{R}_{pde}(P_{\theta},Q_{\theta})+\omega_{1}\mathcal{R}_{b_{1}}(P_{\theta})+\omega_{2}\mathcal{R}_{b_{2}}(P_{\theta},Q_{\theta})$
(9)

where

$\mathcal{R}_{pde}(P_{\theta},Q_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},\delta_{i},\cdot)\right\|^{2}_{L^{2}([0,L]\times[0,T])},$
$\mathcal{R}_{b_{1}}(P_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|P_{\theta}(s_{i},\delta_{i},(L,\cdot))-p_{i}\right\|_{L^{2}([0,T])}^{2},$

and

$\mathcal{R}_{b_{2}}(P_{\theta},Q_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|W(P_{\theta},Q_{\theta})(s_{i},\delta_{i},(L,\cdot))\right\|_{L^{2}([0,T])}^{2}.$

The complete structure of the BP-DeepONet is displayed in Figure 9.

Figure 9: Structure of the BP-DeepONet

To evaluate the different residual terms $\mathcal{R}_{pde}$,
$\mathcal{R}_{b1}$, and $\mathcal{R}_{b2}$, we need to approximate the $L_{2}$
norm with some discrete quadrature. For the PDE residual term
$\mathcal{R}_{pde}$, we can simply approximate the $i$-th $L_{2}$ residual by

$\displaystyle\left\|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},\delta_{i},\cdot)\right\|^{2}_{L^{2}([0,L]\times[0,T])}$
$\displaystyle\approx$
$\displaystyle\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},\delta_{i},(z_{j},t_{j}))\right|^{2},$

where $\\{(z_{j},t_{j})\\}_{j=1}^{M}$ is a set of quadrature points randomly
sampled from the domain $[0,L]\times[0,T]$. The partial derivatives in the
residual terms can be calculated as

$\frac{\partial}{\partial
x}P_{\theta}(s_{i},\delta_{i},(z_{j},t_{j}))=\mathbf{b_{1}}(s_{i})^{\top}\frac{\partial}{\partial
x}\mathbf{k}(\delta_{i},(z_{j},t_{j})),\quad x=z\text{ or }t,$
$\frac{\partial}{\partial
x}Q_{\theta}(s_{i},\delta_{i},(z_{j},t_{j}))=\mathbf{b_{2}}(s_{i})^{\top}\frac{\partial}{\partial
x}\mathbf{k}(\delta_{i},(z_{j},t_{j})),\quad x=z\text{ or }t.$

However, evaluating this approximation is very time-consuming because the
total number of gradient evaluations in one training epoch is $O(MN)$. It is
necessary to design a more efficient method to approximate the residuals. This
can be achieved by utilizing the periodicity of our BP-DeepONet.

###  4.3 Efficient implementation of training BP-DeepONet

Given any periodic function $f(t)$ with period $\delta$, we have the norm
$\|\cdot\|_{L^{2}([0,T])}$ is equivalent to $\|\cdot\|_{L^{2}([0,\delta])}$:

$\left\lfloor\frac{T}{\delta}\right\rfloor\|f\|^{2}_{L^{2}([0,\delta])}\leq\|f\|^{2}_{L^{2}([0,T])}\leq\left\lceil\frac{T}{\delta}\right\rceil\|f\|^{2}_{L^{2}([0,\delta])},$

where $\lfloor\cdot\rfloor$ and $\lceil\cdot\rceil$ denote the floor function
and ceiling function, respectively. Then, for the $i$-th $L_{2}$ residual term
in $\mathcal{R}_{pde}$, we may replace it with

$\displaystyle\left\|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},\delta_{i},\cdot)\right\|^{2}_{L^{2}([0,L]\times[0,\delta_{i}])}$
$\displaystyle=$
$\displaystyle\int_{0}^{\delta_{i}}\int_{0}^{L}\left|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},\delta_{i},(z,t))\right|^{2}dzdt$
$\displaystyle=$
$\displaystyle\int_{0}^{\delta_{i}}\int_{0}^{L}\left|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},1,(z,t/\delta_{i}))\right|^{2}dzdt.$

After a change of variable $\hat{t}=t/\delta_{i}$, the above integration can
be written as

$\displaystyle\int_{0}^{1}\int_{0}^{L}\left|\left(\frac{\partial}{\partial\hat{t}}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}/\delta_{i}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},1,(z,\hat{t}))\right|^{2}dzd\hat{t}$
$\displaystyle\approx$
$\displaystyle\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial\hat{t}}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}/\delta_{i}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},1,(z_{j},\hat{t}_{j}))\right|^{2}$

where $\\{(z_{j},\hat{t}_{j})\\}_{j=1}^{M}$ is randomly sampled from the
domain $[0,L]\times[0,1]$. We approximate $\mathcal{R}_{pde}$ by

$\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial\hat{t}}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}/\delta_{i}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},1,(z_{j},\hat{t}_{j}))\right|^{2}$
(10)

Then, the partial derivatives in (10) can be calculated as

$\frac{\partial}{\partial
x}P_{\theta}(s_{i},\delta_{i},(z_{j},\hat{t}_{j}))=\mathbf{b_{1}}(s_{i})^{\top}\frac{\partial}{\partial
x}\mathbf{k}(1,(z_{j},\hat{t}_{j})),\quad x=z\text{ or }\hat{t},$
$\frac{\partial}{\partial
x}Q_{\theta}(s_{i},1,(z_{j},\hat{t}_{j}))=\mathbf{b_{2}}(s_{i})^{\top}\frac{\partial}{\partial
x}\mathbf{k}(1,(z_{j},\hat{t}_{j})),\quad x=z\text{ or }\hat{t}.$

For different $i$, we do not need to recompute the gradient of $\mathbf{k}$,
so the number of gradient evaluations in one training epoch is only $O(M)$. We
can use the same idea to approximate $\mathcal{R}_{b2}$ by

$\displaystyle\mathcal{R}_{b2}(P_{\theta},Q_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(Q_{\theta}(1+R_{1}/R_{2})+CR_{1}\frac{\partial}{\partial\hat{t}}Q_{\theta}-P_{\theta}/R_{2}-C\frac{\partial}{\partial\hat{t}}P_{\theta}\right)(s_{i},1,(L,\hat{t}_{j})).\right|^{2}$
(11)

For $\mathcal{R}_{b1}$, we simply approximate it by

$\displaystyle\mathcal{R}_{b1}(P_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|P_{\theta}(L,\hat{t}_{j})-\bar{p}_{i}(\hat{t}_{j}\delta_{i})\right|^{2}$
(12)

where $\\{\hat{t}_{j}\\}_{j=1}^{M}$ is randomly sampled from the domain
$[0,1]$ and $\bar{p}_{i}$ is the average of $p_{i}$ over all periods contained
in $[0,T]$. Finally, the entire empirical loss function for training the BP-
DeepONet is the combination of three terms (10), (11), and (12):

$\displaystyle\hat{\mathcal{L}}(P_{\theta},Q_{\theta})=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial\hat{t}}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}/\delta_{i}+H(P_{\theta},Q_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta})\right)(s_{i},1,(z_{j},\hat{t}_{j}))\right|^{2}$
$\displaystyle+\frac{\omega_{1}}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|P_{\theta}(s_{i},1,(L,\hat{t}_{j}))-\bar{p}_{i}(\hat{t}_{j}\delta_{i})\right|^{2}$
$\displaystyle+\frac{\omega_{2}}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(Q_{\theta}(1+R_{1}/R_{2})+CR_{1}\frac{\partial}{\partial\hat{t}}Q_{\theta}-P_{\theta}/R_{2}-C\frac{\partial}{\partial\hat{t}}P_{\theta}\right)(s_{i},1,(L,\hat{t}_{j}))\right|^{2}.$

###  4.4 Meta BP-DeepONet for sample-dependent hyper-parameters estimation

In the Navier-Stokes equation, many hyper-parameters depend on the physical
properties of arterial systems. In the previous BP-DeepONet, we set all these
hyper-parameters to be subject-independent, i.e., we use the same set of
hyper-parameters for all samples. However, in real situations, these
parameters differ from person to person, so it is unrealistic to set them all
as constants. To handle this issue, we generalize the model (5) to be sample-
dependent. For each training sample $i$ from $1$ to $N$, we denote the set of
all hyper-parameters for the $i$-th sample as

$\gamma_{i}:=\\{\beta_{i},A_{0,i},\rho_{i},K_{r,i},P_{ext,i},R_{1,i},R_{2,i},C_{i}\\},$

and we assume the solutions $P_{i},Q_{i}$ for the $i$-th sample satisfy the
following PDE:

$\frac{\partial U_{i}}{\partial t}+H(U_{i},\gamma_{i})\frac{\partial
U_{i}}{\partial z}=B(U_{i},\gamma_{i}),\quad(z,t)\in[0,L]\times[0,T]$ (13)

where

$U_{i}=\begin{pmatrix}P_{i}\\\ Q_{i}\end{pmatrix},$
$H(U_{i},\gamma_{i})=\begin{pmatrix}0&\frac{\beta_{i}}{2a(P_{i},\gamma_{i})A_{0,i}}\\\
\frac{a(P_{i},\gamma_{i})^{2}}{\rho_{i}}-\frac{2A_{0,i}Q_{i}^{2}}{\beta_{i}a^{3}(P_{i},\gamma_{i})}&\frac{2Q_{i}}{a(P_{i},\gamma_{i})^{2}}\end{pmatrix},$
$B(U_{i},\gamma_{i})=\begin{pmatrix}0\\\
-K_{r,i}\frac{Q_{i}}{a(P_{i},\gamma_{i})^{2}}\end{pmatrix},$

and

$a(P_{i},\gamma_{i})=\frac{A_{0,i}}{\beta_{i}}(P_{i}-P_{ext,i})+\sqrt{A_{0,i}}.$

The corresponding boundary and periodic conditions are:

$P_{i}(L,t)=p_{i}(t),\quad t\in[0,T],$
$W(P_{i},Q_{i},\gamma_{i}):=Q_{i}(L,t)(1+\frac{R_{1,i}}{R_{2,i}})+C_{i}R_{1,i}\frac{\partial
Q_{i}}{\partial t}(L,t)-\frac{P_{i}(L,t)}{R_{2,i}}-C_{i}\frac{\partial
P_{i}}{\partial t}(L,t)=0,\quad t\in[0,T],$ (14)

and

$U_{i}(z,t+\delta_{i})=U_{i}(z,t),\quad z\in[0,L].$

It is almost impossible to measure $\gamma_{i}$ for each sample, so we adopt
the idea from the field of meta-learning to estimate $\gamma_{i}$. We further
define a hyper-network on top of the branch net. This network takes
$\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ as inputs and outputs the estimation of
$\gamma_{i}$. This hyper-network is optimized together with the BP-DeepONet.
We also impose extra constraints to regularize the output of the hyper-
network. The first constraint is imposed as:

$\int_{0}^{T}P_{i}(z,t)dt=\int_{0}^{T}p_{i}(t)dt,\quad\forall z\in[0,L],$ (15)

which implies that the mean blood pressure does not change during the wave
propagation. This phenomenon was found in clinical observations [54]. The
second constraint is imposed as

$R_{1,i}+R_{2,i}=\frac{\int_{0}^{\delta_{i}}P_{i}(L,t)dt}{\int_{0}^{\delta_{i}}Q_{i}(L,t)dt}.$
(16)

This equation can be derived by integrating the equation (14) over one period
$[0,\delta_{i}]$ and it is used for regularizing the estimation of $R_{1,i}$
and $R_{2,i}$. If we have other information about the hemodynamics, like the
pulse wave velocity, we may introduce them as constraints to regularize the
hyper-networks. Then, the training loss function of the meta BP-DeepONet is

$\displaystyle\mathcal{L}_{meta}(P_{\theta},Q_{\theta},\gamma_{\theta})=$
$\displaystyle\mathcal{R}_{pde}(P_{\theta},Q_{\theta},\gamma_{\theta})$
$\displaystyle+\omega_{1}\mathcal{R}_{b_{1}}(P_{\theta})+\omega_{2}\mathcal{R}_{b_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})+\omega_{3}\mathcal{R}_{c_{1}}(P_{\theta},Q_{\theta})+\omega_{4}\mathcal{R}_{c_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})$
(17)

where $\Gamma_{\theta}$ is the set of estimated hyper-parameters given by the
hyper-network:

$\gamma_{\theta}(s_{i})=\\{\beta_{\theta}(s_{i}),A_{0,\theta}(s_{i}),\rho_{\theta}(s_{i}),K_{r,\theta}(s_{i}),P_{ext,\theta}(s_{i}),R_{1,\theta}(s_{i}),R_{2,\theta}(s_{i}),C_{\theta}(s_{i})\\},\quad
i=1,\dots,N,$

and the residual terms in (17) are defined as

$\mathcal{R}_{pde}(P_{\theta},Q_{\theta},\gamma_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta},\gamma_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta},\gamma_{\theta})\right)(s_{i},\delta_{i},\cdot)\right\|^{2}_{L^{2}([0,L]\times[0,T])},$
$\mathcal{R}_{b_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|W(P_{\theta},Q_{\theta},\gamma_{\theta})(s_{i},\delta_{i},(L,\cdot))\right\|_{L^{2}([0,T])}^{2},$
$\mathcal{R}_{c_{1}}(P_{\theta},Q_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\int_{0}^{T}P_{\theta}(s_{i},\delta_{i},(\cdot,t))dt-\int_{0}^{T}p_{i}(t)dt\right\|_{L^{2}([0,L])}^{2},$
$\mathcal{R}_{c_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left(R_{1,\theta}(s_{i})+R_{2,\theta}(s_{i})-\frac{\int_{0}^{\delta_{i}}P_{\theta}(s_{i},\delta_{i},(L,t))dt}{\int_{0}^{\delta_{i}}Q_{\theta}(s_{i},\delta_{i},(L,t))dt}\right)^{2},$

and $\mathcal{R}_{b_{1}}(P_{\theta})$ is defined the same as (9). To
approximate these residual terms, we use the same change of variable technique
as section 4.2:

$\mathcal{R}_{pde}(P_{\theta},Q_{\theta},\gamma_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta},\gamma_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta},\gamma_{\theta})\right)(s_{i},1,(z_{j},\hat{t}_{j}))\right|^{2},$
$\mathcal{R}_{b_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(Q_{\theta}(1+R_{1,\theta}/R_{2,\theta})+C_{\theta}R_{1,\theta}\frac{\partial}{\partial\hat{t}}Q_{\theta}-P_{\theta}/R_{2,\theta}-C_{\theta}\frac{\partial}{\partial\hat{t}}P_{\theta}\right)(s_{i},1,(L,\hat{t}_{j}))\right|^{2},$
$\mathcal{R}_{c_{1}}(P_{\theta},Q_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{j_{1}=1}^{M^{\prime}}\left|\sum_{j_{2}=1}^{M}(P_{\theta}(s_{i},1,(\hat{z}_{j_{1}},\hat{t}_{j_{2}}))-p_{i}(\hat{t}_{j_{2}}))\right|^{2},$
$\mathcal{R}_{c_{2}}(P_{\theta},Q_{\theta},\gamma_{\theta})\approx\frac{1}{N}\sum_{i=1}^{N}\left|R_{1,\theta}(s_{i})+R_{2,\theta}(s_{i})-\frac{\sum_{j=1}^{M}P_{\theta}(s_{i},\delta_{i},(L,\hat{t}_{j}))}{\sum_{j=1}^{M}Q_{\theta}(s_{i},\delta_{i},(L,\hat{t}_{j}))}\right|^{2},$

where $\\{(z_{j},\hat{t}_{j})\\}_{j=1}^{M}$ is a set of quadrature points
sampled from $[0,L]\times[0,1]$ and $\\{\hat{z}_{j}\\}_{j=1}^{M^{\prime}}$ is
a set of quadrature points sampled from $[0,L]$, and
$\mathcal{R}_{b_{1}}(P_{\theta})$ is approximated in the same way as (12).
Then, the empirical training loss for the meta BP-DeepONet is given as

$\displaystyle\hat{\mathcal{L}}_{meta}(P_{\theta},Q_{\theta},\gamma_{\theta})=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(\frac{\partial}{\partial
t}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}+H(P_{\theta},Q_{\theta},\gamma_{\theta})\frac{\partial}{\partial
z}\begin{pmatrix}P_{\theta}\\\
Q_{\theta}\end{pmatrix}-B(P_{\theta},Q_{\theta},\gamma_{\theta})\right)(s_{i},1,(z_{j},\hat{t}_{j}))\right|^{2}$
$\displaystyle+\frac{\omega_{1}}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|P_{\theta}(s_{i},1,(L,\hat{t}_{j}))-\bar{p}_{i}(\hat{t}_{j}\delta_{i})\right|^{2}$
$\displaystyle+\frac{\omega_{2}}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left|\left(Q_{\theta}(1+R_{1,\theta}/R_{2,\theta})+C_{\theta}R_{1,\theta}\frac{\partial}{\partial\hat{t}}Q_{\theta}-P_{\theta}/R_{2,\theta}-C_{\theta}\frac{\partial}{\partial\hat{t}}P_{\theta}\right)(s_{i},1,(L,\hat{t}_{j}))\right|^{2}$
$\displaystyle+\frac{\omega_{3}}{N}\sum_{i=1}^{N}\sum_{j_{1}=1}^{M^{\prime}}\left|\sum_{j_{2}=1}^{M}(P_{\theta}(s_{i},1,(\hat{z}_{j_{1}},\hat{t}_{j_{2}}))-p_{i}(\hat{t}_{j_{2}}))\right|^{2}$
$\displaystyle+\frac{\omega_{4}}{N}\sum_{i=1}^{N}\left|R_{1,\theta}(s_{i})+R_{2,\theta}(s_{i})-\frac{\sum_{j=1}^{M}P_{\theta}(s_{i},\delta_{i},(L,\hat{t}_{j}))}{\sum_{j=1}^{M}Q_{\theta}(s_{i},\delta_{i},(L,\hat{t}_{j}))}\right|^{2}.$

The overall structures of the meta BP-DeepONet are given in Figure 10 where
the hyper-network is implemented as a simple FCN.

Figure 10: Strucutres of the meta BP-DeepONet

######  Remark.

During the training of the meta BP-DeepONet, we only require the measurement
of blood pressure at $z=L$. Then the whole solution, all boundary conditions,
initial conditions, and hyper-parameters can be automatically learned from
this measurement. No extra measurements and parameters are needed.

##  5 Numerical experiments of the PINN on simulated data

In this section, we would like to validate the proposed physics-informed
training method by solving one PDE instance. We first generate a solution to
the Navier-Stokes equation (4) using a MacCormack scheme with a prescribed
initial condition and inlet boundary condition. Then, we solve this equation
with the boundary conditions (7) and (6) and the time-periodic condition (8)
using a simple PINN. We can see in the numerical results that the PINN
solution is very close to the simulated solution.

We consider the 1D Navier-stokes equation (1) and (2) with the initial
condition:

$A(x,0)=A_{0},\quad Q(x,0)=0,$

the inflow boundary condition:

$Q(0,t)=q_{0}(t),$

and the three-element Windkessel boundary condition

$Q(L,t)(1+\frac{R_{1}}{R_{2}})+CR_{1}\frac{\partial Q}{\partial
t}(L,t)=\frac{P(L,t)}{R_{2}}+C\frac{\partial P}{\partial t}(L,t).$

Such an initial boundary value problem can be efficiently solved by many
classical methods like the MacCormack scheme[26, 27], Taylor-Galerkin scheme
[28, 29], MUSCL [30, 31], and the local discontinuous Galerkin scheme [32, 33,
34]. We use a numerical example that is similar to the baseline aorta case
provided in [21], where the hyper-parameters are chosen as $L=25\text{cm}$,
$A_{0}=4.52\text{cm}^{2}$, $\beta=1134.37\text{kg/s}^{2}$,
$\rho=1060\text{kg/m}^{3}$, $P_{ext}=9\text{kPa}$, $R_{1}=1.17\times
10^{7}\text{Pa s m}^{-3}$, $R_{1}=1.12\times 10^{8}\text{Pa s m}^{-3}$,
$C=1.01\times 10^{-8}\text{m}^{3}\text{Pa}^{-1}$, $v=0.01\times
10^{6}\text{m}^{2}/\text{s}^{2}$, and the inflow boundary conditions
$q_{0}(t)$ is given in Figure 11. Then, we solve this problem using the
MacCormack scheme.

Figure 11: The inflow boundary condition

After having a referenced numerical solution $\hat{A}$ and $\hat{Q}$, we
define a residual network $N_{\theta}(z,t)=(P_{\theta}(z,t),Q_{\theta}(z,t))$
shown in Figure 4 (b) with five residual layers. The time-periodic condition
is incorporated into the network using the same feature expansion technique
described before. We train this network by minimizing the following loss
function:

$\displaystyle L_{PINN}$ $\displaystyle(\theta):=\left\|\frac{\partial
A_{\theta}}{\partial t}+\frac{\partial Q_{\theta}}{\partial
z}\right\|_{L^{2}([0,L]\times[0,T])}^{2}+\left\|\frac{\partial
Q_{\theta}}{\partial t}+\frac{\partial}{\partial
z}(\frac{Q_{\theta}^{2}}{A_{\theta}})+\frac{A_{\theta}}{\rho}\frac{\partial
P_{\theta}}{\partial
z}+K_{r}\frac{Q_{\theta}}{A_{\theta}}\right\|_{L^{2}([0,L]\times[0,T])}^{2}$
$\displaystyle+100\left(\|P_{\theta}(L,\cdot)-\hat{P}(L,\cdot)\|_{L^{2}([0,T])}^{2}+\left\|Q_{\theta}(L,\cdot)(1+\frac{R_{1}}{R_{2}})+CR_{1}\frac{\partial
Q_{\theta}}{\partial
t}(L,\cdot)-\frac{P_{\theta}(L,\cdot)}{R_{2}}-C\frac{\partial
P_{\theta}}{\partial t}(L,\cdot)\right\|_{L^{2}([0,T])}^{2}\right)$

where
$P_{\theta}(z,t)=P_{ext}+\frac{\beta}{A_{0}}(\sqrt{A_{\theta}(z,t)}-\sqrt{A_{0}})$.
We use a standard Adam optimizer with a learning rate of 1e-4, and the total
number of training epochs is 50,000. We randomly sample 1,000 quadrature
points in the domain in each epoch to approximate the $L^{2}$ integration. We
can see the comparison of the PINN solution and the numerical solution by the
MacCormack scheme in Figure 12. The PINN solution solved using the loss
function $L_{PINN}$ is very close to the solution of the MacCormack scheme.

Figure 12: Comparison of the PINN solution and the finite difference solution
(MacCormack).

We also calculate the relative error:

$\left\|\frac{\hat{P}-P_{\theta}}{\hat{P}}\right\|_{1}\approx 0.30\%$

and

$\frac{\|\hat{Q}-Q_{\theta}\|_{1}}{\|\hat{Q}\|_{\infty}}\approx 0.67\%.$

These two metrics are also used in [21] to measure the error of their
numerical solutions for the Navier-Stokes equation. From Figure 12, we can
observe that the error of the PINN solution is larger near the dicrotic notch
where the solution is near singular. This case is very common when using the
standard PINN method to solve PDEs, and it can be improved by using some
adaptive sampling scheme [55]. We may consider this direction in our future
work.

##  6 Numerical experiments of the PI-DeepONet methods on clinical data

###  6.1 Dataset description

We conducted our experiments on the MIMIC dataset [56], the largest and most
commonly used dataset for cuffless blood pressure estimation. We use the pre-
processed version of the MIMIC data from [57], which contains recordings of
12,000 subjects after simple pre-processing. Each recording consists of
continuous ECG, PPG, and ABP waveform data. The length of each recording is
within 10 minutes. All these data are collected from ICU patients, and the
invasive method collects the ABP waveform data. Furthermore, we also apply
some extra pre-processing steps to these data:

  1. 1.

Use a bandpass filter to remove noise from PPG and ABP signals.

  2. 2.

Align the PPG and ABP signals.

  3. 3.

Randomly select 75% of subjects for training and 25% for testing.

  4. 4.

20% of the training set is selected as the validation set.

  5. 5.

The first 15% of data in each testing subject is used for calibration.

  6. 6.

Split each recording into many 4-second segments without overlapping.

  7. 7.

Remove those segments that are highly non-periodic.

  8. 8.

Remove subjects whose number of valid segments is less than 10.

  9. 9.

Apply some peak detection algorithms to detect all peaks in each PPG segment,
then calculate the average length of intervals between two consecutive peaks.

After pre-processing, the total number of training and testing samples is
about 70,000. Step 5 is also called individual fine-tuning [58], and it is a
necessary process to calibrate and improve the blood pressure measurement. The
segments for calibration will not used when calculating the testing error. One
segment of the recording is shown in Figure 13.

Figure 13: One segment of ECG, PPG, and ABP signals.

###  6.2 Training procedures

For the training of the BP-DeepONet and the Meta BP-DeepONet, We apply the
standard Adam optimizer with a scheduled learning rate to train the network.
The initial learning rate is set to be 1e-4, and decays by a factor of 0.9
every 5 epochs. After 200 epochs, we restart the learning rate with 1e-4 and
continue to train with the same learning rate schedule. After several rounds
of learning rate cycles, we choose the best model based on the performance
evaluated on a validation set.

When training the BP-DeepONet, we need first to determine all the hyper-
parameters in the NS model, so we set hyper-parameters based on the reference
values of the radial artery described in [59] and estimate the Windkessel
parameters using the method provided in [21]: $A_{0}=8\times 10^{-6}\ m^{2}$,
$\beta=529\ kg/s^{2}$, $L=0.25m$, $K_{r}=8*\pi*0.01\approx 0.251$, $\rho=1.06\
kg/m^{3}$, $P_{ext}=6.67\ kPa$, $R_{1}=12.37\ Pa\cdot s/m^{3}$, $R_{2}=137.5\
Pa\cdot s/m^{3}$, $C=0.01\ m^{3}/Pa$. In the Meta BP-DeepONet, we also use
these values as referenced values for the predicted hyper-parameters. We
restrict the deviation of the predicted hyper-parameters from the referenced
values to be less than or equal to 20 % so that the estimation of the hyper-
parameters will always be reasonable.

During the training of both BP-DeepONet and meta BP-DeepONet, we use a batch
size of $100$. All the training is done on a GPU server with 4 NVIDIA RTX
A6000 GPUs.

###  6.3 Numerical results

After finishing the training of DeepONet, we can obtain the prediction of the
BP waveform at any given location $z$ by evaluating the DeepONet at
$(z,t_{0})$, $(z,t_{1})$, …, $(z,t_{N_{h}})$, where
$\\{t_{j}\\}_{j=0}^{N_{h}}$ is a uniform discretization of the domain $[0,T]$.

We first evaluate the waveform prediction at the outlet side $z=L$. Since we
have the measured BP wave at $z=L$ by the invasive method, we can compute the
prediction error of our method at this location. We evaluate the error using
several metrics: waveform MAE, SBP MAE, DBP MAE, and MBP MAE, where MAE is
short for mean absolute error. Suppose
$P_{i}=(P_{\theta}(s_{i},\delta_{i},(L,t_{0})),P_{\theta}(s_{i},\delta_{i},(L,t_{1})),\dots,P_{\theta}(s_{i},\delta_{i},(L,t_{N_{h}})))\in\mathbb{R}^{N_{h}}$
is the predicted waveform and $p_{i}\in\mathbb{R}^{N_{h}}$ is the ground true
BP waveform. Then, the four metrics are defined as follows:

  * 1.

Waveform MAE: $\frac{1}{N}\frac{1}{N_{h}}\sum_{i=1}^{N}|P_{i}-p_{i}|_{1}$,

  * 2.

SBP MAE: $\frac{1}{N}\sum_{i=1}^{N}|\max(P_{i})-\max(p_{i})|$,

  * 3.

DBP MAE: $\frac{1}{N}\sum_{i=1}^{N}|\min(P_{i})-\min(p_{i})|$,

  * 4.

MBP MAE: $\frac{1}{N}\sum_{i=1}^{N}|\text{mean}(P_{i})-\text{mean}(p_{i})|$,

where $|\cdot|_{1}$ denotes the standard vector $l_{1}$ norm. In Table 1, we
listed the error metrics for both BP-DeepONet and meta BP-DeepONet. We also
compared our results with the IEEE standard for wearable, cuffless blood
pressure measuring devices [60]. The suggested grading is also shown in Table
1. We can see the meta BP-DeepONet achieves better accuracy on the MAE of SBP,
DBP, and MBP, while the BP-DeepONet achieves better accuracy on the waveform
MAE. According to the IEEE standard, the meta BP-DeepONet is graded as A-level
in all metrics, while the BP-DeepONet is graded as A-level except on the SBP
measurement.

|  Methods |  Waveform |  SBP |  DBP |  MBP
---|---|---|---|---|---
MAE (mmHg) |  BP-DeepONet |  4.406 |  5.057 |  2.581 |  2.881
Meta BP-DeepONet |  4.555 |  4.934 |  2.574 |  2.871
Grading |  BP-DeepONet |  A |  B |  A |  A
Meta BP-DeepONet |  A |  A |  A |  A
Table 1: The prediction error of two proposed methods and their grading
according to the IEEE standard.

Since our methods can predict the blood pressure $P(z,t)$ and blood flow rate
$Q(z,t)$ in the whole domain $(z,t)\in[0,L]\times[0,T]$, we also want to
evaluate the errors in the interior of this domain. We choose two different
metrics to evaluate the error in $(z,t)\in[0,L]\times[0,T]$. The first metric
is the residual errors $\mathcal{R}_{pde}$. The residual error is usually used
to estimate the total error when the true solution is unknown. For the BP-
DeepONet, the residual error is defined as:

$\displaystyle
RE_{1}(P_{\theta},Q_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial
A_{\theta}}{\partial t}+\frac{\partial Q_{\theta}}{\partial
z}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{2}([0,L]\times[0,T])}$
$\displaystyle
RE_{2}(P_{\theta},Q_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial
Q_{\theta}}{\partial t}+\frac{\partial}{\partial
z}(\frac{Q_{\theta}^{2}}{A_{\theta}})+\frac{A_{\theta}}{\rho}\frac{\partial
P_{\theta}}{\partial
z}+K_{r}\frac{Q_{\theta}}{A_{\theta}}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{2}([0,L]\times[0,T])}$

where
$A_{\theta}=\left(\frac{A_{0}}{\beta}(P_{\theta}-P_{ext})+\sqrt{A_{0}}\right)^{2}$.
For the meta BP-DeepONet, the residual error is defined as:

$\displaystyle
RE_{1}(P_{\theta},Q_{\theta},\gamma_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial
A_{\theta,i}}{\partial t}+\frac{\partial Q_{\theta}}{\partial
z}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{2}([0,L]\times[0,T])}$
$\displaystyle
RE_{2}(P_{\theta},Q_{\theta},\gamma_{\theta})=\frac{1}{N}\sum_{i=1}^{N}\left\|\left(\frac{\partial
Q_{\theta}}{\partial t}+\frac{\partial}{\partial
z}(\frac{Q_{\theta}^{2}}{A_{\theta,i}})+\frac{A_{\theta,i}}{\rho_{i}}\frac{\partial
P_{\theta}}{\partial
z}+K_{r_{i}}\frac{Q_{\theta}}{A_{\theta,i}}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{2}([0,L]\times[0,T])}$

where
$A_{\theta,i}=\left(\frac{A_{0,i}}{\beta}(P_{\theta}-P_{ext,i})+\sqrt{A_{0,i}}\right)^{2}$.
In [43], the authors show that the error of solutions to the Navier-Stokes by
neural networks can be bounded by the residual errors. The residual errors of
the two methods are shown in Table 2.

|  $RE_{1}(\sqrt{L})$ |  $RE_{2}(\sqrt{dm*L/s})$
---|---|---
BP-DeepONet |  1.11$\times$ 1e-4 |  2.65$\times$ 1e-2
Meta BP-DeepONet |  0.22$\times$ 1e-4 |  1.14$\times$1e-2
Table 2: The residual error of two proposed methods.

The second metric we used to evaluate the error is the relative error compared
to a fine-tuned BP-DeepONet (or meta BP-DeepONet) on each test sample. Though
we do not know the ground true solution, we have shown in Section 5 that a
well-trained PINN solution is reasonably accurate compared to the simulated
true solution. Therefore, we can fine-tune the network on each test sample,
i.e., minimize the PDE residual loss of this particular sample only, and use
this fine-tuned network as an approximation to the true solution. During the
fine-tuning stage, we fix the parameters in the branch net and only train the
parameters in the trunk net so that predictions of hyper-parameters will not
change. The error metric is defined as:

$\sum_{i=1}^{N}\frac{\left\|\left(\hat{P}_{tuned}-P_{\theta}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{1}([0,L]\times[0,T])}}{\left\|\hat{P}_{tuned}(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{1}([0,L]\times[0,T])}}$

and

$\sum_{i=1}^{N}\frac{\left\|\left(\hat{Q}_{tuned}-Q_{\theta}\right)(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{1}([0,L]\times[0,T])}}{\left\|\hat{Q}_{tuned}(s_{i},\delta_{i},(\cdot,\cdot))\right\|_{L^{\infty}([0,L]\times[0,T])}}$

where $\hat{P}_{tuned}$ and $\hat{Q}_{tuned}$ denote the solution after fine-
tuning. The relative errors are given in Table 3. We can see the BP-DeepONet
can predict the PDE solution better than the Meta BP-DeepONet. This result is
not surprising because the BP-DeepONet assumes all samples satisfy the same
PDE, which greatly simplifies the problem.

|  error of $P_{\theta}$ |  error of $Q_{\theta}$
---|---|---
BP-DeepONet |  4.10$\times$ 1e-2 |  5.35$\times$ 1e-2
Meta BP-DeepONet |  4.18$\times$ 1e-2 |  1.03$\times$1e-1
Table 3: The relative error of two proposed methods.

We also visualize the predicted waveforms for one test sample in Figure 14 and
Figure 15. We can observe that the predictions of both the BP-DeepONet and the
meta BP-DeepONet are very close to the ground truth data at the outlet side
$z=L$ (Figure 14 (c)). We can also observe the pulse pressure amplification
with the difference between the maximum and minimum of the waveforms increases
when $z$ increases for both methods in Figure 15. Because the Meta BP-DeepONet
uses different hyper-parameters for each sample, its prediction differs from
BP-DeepONet in $z\in(0,1)$, especially at $z=0$.

Figure 14: (a) and (b) show the input physiological signals. (c) shows the
prediction and ground truth data at $z=L$.

Figure 15: Predictions of ABP waveforms at different locations $z$.

##  7 Conclusions and Future Works

In this work, we propose two methods based on the physics-informed DeepONet to
predict the blood pressure waveform: the BP-DeepONet and the meta BP-DeepONet.
We incorporate the Navier-Stokes equation with time-periodic conditions and
Windkessel boundary conditions into the physics-informed training procedure
and provide efficient implementation algorithms. The proposed methods are the
first ones that can predict waveforms at different locations. During our
numerical experiments, the proposed methods can predict blood pressure
waveforms at reasonable accuracy and preserve some physical properties
consistent with the clinical observations.

This work only considers very standard DNN models for trunk and branch nets.
We may try implementing different DNN architectures in the future, like the
transformer, which is very powerful in processing time-series data. Besides,
the inputs to the DeepONet are just the pre-processed physiological signals.
We may also try including some morphological features in the inputs,
effectively improving the prediction accuracy.

We know the arterial system of humans is very complicated. As the first work
in this direction, we only consider a simple cylindrical segment of an artery.
Still, we may generalize our methods to more complex geometry afterward—for
example, a simple arterial tree with bifurcations.

In the future, we will also consider generating a simulated dataset to
validate better the accuracy of both BP-DeepONet and meta BP-DeepONet methods.

## Acknowledgement

The work of Lingfeng Li is supported by the InnoHK project at Hong Kong Centre
for Cerebro-cardiovascular Health Engineering (COCHE). The work of Xue-Cheng
Tai is partially supported by HKRGC-NSFC Grant N-CityU214/19, HKRGC CRF. Grant
C1013-21GF and NORCE Kompetanseoppbygging program. The work of Raymond Chan is
partially supported by HKRGC GRF grants CityU1101120, CityU11309922, CRF grant
C1013-21GF and HKITF MHKJFS Grant MHP/054/22. We want to thank Prof Yuan-Ting
Zhang and Dr. Alireza Keramat for the helpful discussion during this project.

## References

  * [1] World Health Organization, Cardiovascular diseases (cvds), https://www.who.int/news-room/fact-sheets/detail/cardiovascular-diseases-(cvds).
  * [2] F. D. Fuchs, P. K. Whelton, High blood pressure and cardiovascular disease, Hypertension 75 (2) (2020) 285–292.
  * [3] B. Matthieu, N.-G. Karine, C. Vincent, V. Alain, C. J. François, R. Philippe, S. François, Cardiac output measurement in patients undergoing liver transplantation: pulmonary artery catheter versus uncalibrated arterial pressure waveform analysis, Anesthesia & Analgesia 106 (5) (2008) 1480–1486.
  * [4] R. H. Thiele, M. E. Durieux, Arterial waveform analysis for the anesthesiologist: past, present, and future concepts, Anesthesia & Analgesia 113 (4) (2011) 766–776.
  * [5] Y. Li, H. Gu, H. Fok, J. Alastruey, P. Chowienczyk, Forward and backward pressure waveform morphology in hypertension, Hypertension 69 (2) (2017) 375–381.
  * [6] D. A. Hullender, O. R. Brown, Simulations of blood pressure and identification of atrial fibrillation and arterial stiffness using an extended kalman filter with oscillometric pulsation measurements, Computer Methods and Programs in Biomedicine 198 (2021) 105768.
  * [7] P. Su, X.-R. Ding, Y.-T. Zhang, J. Liu, F. Miao, N. Zhao, Long-term blood pressure prediction with deep recurrent neural networks, in: 2018 IEEE EMBS International Conference on Biomedical & Health Informatics (BHI), IEEE, 2018, pp. 323–328.
  * [8] H. Eom, D. Lee, S. Han, Y. S. Hariyani, Y. Lim, I. Sohn, K. Park, C. Park, End-to-end deep learning architecture for continuous blood pressure estimation using attention mechanism, Sensors 20 (8) (2020) 2338.
  * [9] N. Ibtehaz, S. Mahmud, M. E. Chowdhury, A. Khandakar, M. Salman Khan, M. A. Ayari, A. M. Tahir, M. S. Rahman, Ppg2abp: Translating photoplethysmogram (ppg) signals to arterial blood pressure (abp) waveforms, Bioengineering 9 (11) (2022) 692.
  * [10] C. Ma, P. Zhang, F. Song, Y. Sun, G. Fan, T. Zhang, Y. Feng, G. Zhang, Kd-informer: Cuff-less continuous blood pressure waveform estimation approach based on single photoplethysmography, IEEE Journal of Biomedical and Health Informatics (2022).
  * [11] K. Lakhal, V. Robert-Edan, Invasive monitoring of blood pressure: a radiant future for brachial artery as an alternative to radial artery catheterisation?, Journal of Thoracic Disease 9 (12) (2017) 4812.
  * [12] S.-P. E, Continuous non-invasive arterial pressure shows high accuracy in comparison to invasive intra-arterial blood pressure measurement, 2008.
URL https://api.semanticscholar.org/CorpusID:30046018

  * [13] M. K. Armstrong, M. G. Schultz, D. S. Picone, J. A. Black, N. Dwyer, P. Roberts-Thomson, J. E. Sharman, Brachial and radial systolic blood pressure are not the same: evidence to support the popeye phenomenon, Hypertension 73 (5) (2019) 1036–1041.
  * [14] J. Kyriazis, J. Glotsos, L. Bilirakis, N. Smirnioudis, The (dp/dt) max derived from arterial pulse waveforms: prospective applications in the haemodialysis setting, Nephrology Dialysis Transplantation 16 (5) (2001) 1087–1088.
  * [15] J. Sharman, A. Qasem, L. Hanekom, D. Gill, R. Lim, T. Marwick, Radial pressure waveform dp/dt max is a poor indicator of left ventricular systolic function, European Journal of Clinical Investigation 37 (4) (2007) 276–281.
  * [16] S. Wang, H. Wang, P. Perdikaris, Learning the solution operator of parametric partial differential equations with physics-informed deeponets, Science Advances 7 (40) (2021) eabi8605.
  * [17] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature Machine Intelligence 3 (3) (2021) 218–229.
  * [18] 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.
  * [19] L. Formaggia, D. Lamponi, A. Quarteroni, One-dimensional models for blood flow in arteries, Journal of Engineering Mathematics 47 (2003) 251–276.
  * [20] A. Quarteroni, L. Formaggia, Mathematical modelling and numerical simulation of the cardiovascular system, Handbook of Numerical Analysis 12 (2004) 3–127.
  * [21] N. Xiao, J. Alastruey, C. Alberto Figueroa, A systematic comparison between 1-d and 3-d hemodynamics in compliant arterial models, International Journal for Numerical Methods in Biomedical Engineering 30 (2) (2014) 204–231.
  * [22] J. Alastruey, A. W. Khir, K. S. Matthys, P. Segers, S. J. Sherwin, P. R. Verdonck, K. H. Parker, J. Peiró, Pulse wave propagation in a model human arterial network: assessment of 1-d visco-elastic simulations against in vitro measurements, Journal of Biomechanics 44 (12) (2011) 2250–2258.
  * [23] M. Saito, Y. Ikenaga, M. Matsukawa, Y. Watanabe, T. Asada, P.-Y. Lagree, One-dimensional model for propagation of a pressure wave in a model of the human arterial network: Comparison of theoretical and experimental results, Journal of Biomechanical Engineering 133 (12) (2011).
  * [24] M. S. Olufsen, C. S. Peskin, W. Y. Kim, E. M. Pedersen, A. Nadim, J. Larsen, Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions, Annals of Biomedical Engineering 28 (2000) 1281–1299.
  * [25] B. N. Steele, J. Wan, J. P. Ku, T. J. Hughes, C. A. Taylor, In vivo validation of a one-dimensional finite-element method for predicting blood flow in cardiovascular bypass grafts, IEEE Transactions on Biomedical Engineering 50 (6) (2003) 649–656.
  * [26] D. Elad, D. Katz, E. Kimmel, S. Einav, Numerical schemes for unsteady fluid flow through collapsible tubes, Journal of Biomedical Engineering 13 (1) (1991) 10–18.
  * [27] J.-M. Fullana, S. Zaleski, A branched one-dimensional model of vessel networks, Journal of Fluid Mechanics 621 (2009) 183–204.
  * [28] J. Wan, B. Steele, S. A. Spicer, S. Strohband, G. R. Feijo´ o, T. J. Hughes, C. A. Taylor, A one-dimensional finite element method for simulation-based medical planning for cardiovascular disease, Computer Methods in Biomechanics & Biomedical Engineering 5 (3) (2002) 195–206.
  * [29] S. J. Sherwin, L. Formaggia, J. Peiro, V. Franke, Computational modelling of 1d blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system, International Journal for Numerical Methods in Fluids 43 (6-7) (2003) 673–700.
  * [30] N. Cavallini, V. Caleffi, V. Coscia, Finite volume and weno scheme in one-dimensional vascular system modelling, Computers & Mathematics with Applications 56 (9) (2008) 2382–2397.
  * [31] O. Delestre, P.-Y. Lagrée, A ‘well-balanced’finite volume scheme for blood flow simulation, International Journal for Numerical Methods in Fluids 72 (2) (2013) 177–205.
  * [32] K. S. Matthys, J. Alastruey, J. Peiró, A. W. Khir, P. Segers, P. R. Verdonck, K. H. Parker, S. J. Sherwin, Pulse wave propagation in a model human arterial network: assessment of 1-d numerical simulations against in vitro measurements, Journal of Biomechanics 40 (15) (2007) 3476–3486.
  * [33] E. Marchandise, M. Willemet, V. Lacroix, A numerical hemodynamic tool for predictive vascular surgery, Medical Engineering & Physics 31 (1) (2009) 131–144.
  * [34] J. Mynard, P. Nithiarasu, A 1d arterial blood flow model incorporating ventricular pressure, aortic valve and regional coronary flow using the locally conservative galerkin (lcg) method, Communications in Numerical Methods in Engineering 24 (5) (2008) 367–417.
  * [35] X. Wang, 1d modeling of blood flow in networks: Numerical computing and applications, Ph.D. thesis, Paris 6 (2014).
  * [36] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895 (2020).
  * [37] L. Lu, X. Meng, Z. Mao, G. E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM Review 63 (1) (2021) 208–228.
  * [38] B. Yu, et al., The deep ritz method: a deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics 6 (1) (2018) 1–12.
  * [39] L. Lyu, Z. Zhang, M. Chen, J. Chen, MIM: A deep mixed residual method for solving high-order partial differential equations, Journal of Computational Physics 452 (2022) 110930.
  * [40] J. Yang, Q. Zhu, A local deep learning method for solving high order partial differential equations, Numerical Mathematics-Theory Methods and Applications (2021).
  * [41] Y. Lu, J. Lu, M. Wang, A priori generalization analysis of the deep ritz method for solving high dimensional elliptic partial differential equations, in: Conference on Learning Theory, PMLR, 2021, pp. 3196–3241.
  * [42] J. Müller, M. Zeinhofer, Error estimates for the variational training of neural networks with boundary penalty, arXiv preprint arXiv:2103.01007 (2021).
  * [43] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes, IMA Journal of Numerical Analysis 42 (2) (2022) 981–1022.
  * [44] L. Li, X.-C. Tai, J. Yang, Generalization error analysis of neural networks with gradient based regularization, Communications in Computational Physics 32 (4) (2022) 1007–1038.
  * [45] L. Li, X.-c. Tai, J. Yang, Q. Zhu, Priori error estimate of deep mixed residual method for elliptic pdes, arXiv preprint arXiv:2206.07474 (2022).
  * [46] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 770–778.
  * [47] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, G. E. Karniadakis, A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114778.
  * [48] Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, A. Anandkumar, Physics-informed neural operator for learning partial differential equations, arXiv preprint arXiv:2111.03794 (2021).
  * [49] S. Goswami, M. Yin, Y. Yu, G. E. Karniadakis, A physics-informed variational deeponet for predicting crack path in quasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114587.
  * [50] N. Westerhof, J.-W. Lankhaar, B. E. Westerhof, The arterial windkessel, Medical & Biological Engineering & Computing 47 (2) (2009) 131–141.
  * [51] J. Alastruey, K. H. Parker, S. J. Sherwin, et al., Arterial pulse wave haemodynamics, in: 11th International Conference on Pressure Surges, Vol. 30, Virtual PiE Led t/a BHR Group Lisbon, Portugal, 2012, pp. 401–443.
  * [52] A. Paviglianiti, V. Randazzo, S. Villata, G. Cirrincione, E. Pasero, A comparison of deep learning techniques for arterial blood pressure prediction, Cognitive Computation 14 (5) (2022) 1689–1710.
  * [53] L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, S. G. Johnson, Physics-informed neural networks with hard constraints for inverse design, SIAM Journal on Scientific Computing 43 (6) (2021) B1105–B1132.
  * [54] A. P. Avolio, L. M. Van Bortel, P. Boutouyrie, J. R. Cockcroft, C. M. McEniery, A. D. Protogerou, M. J. Roman, M. E. Safar, P. Segers, H. Smulyan, Role of pulse pressure amplification in arterial hypertension: experts’ opinion and review of the data, Hypertension 54 (2) (2009) 375–383.
  * [55] Z. Gao, L. Yan, T. Zhou, Failure-informed adaptive sampling for pinns, SIAM Journal on Scientific Computing 45 (4) (2023) A1971–A1994.
  * [56] A. E. Johnson, T. J. Pollard, L. Shen, L.-w. H. Lehman, M. Feng, M. Ghassemi, B. Moody, P. Szolovits, L. Anthony Celi, R. G. Mark, Mimic-iii, a freely accessible critical care database, Scientific Data 3 (1) (2016) 1–9.
  * [57] M. Kachuee, M. Kiani, M. Hoda, M. Shabany, Cuff-Less Blood Pressure Estimation, UCI Machine Learning Repository, DOI: https://doi.org/10.24432/C5B602 (2015).
  * [58] J. Hong, J. Gao, Q. Liu, Y. Zhang, Y. Zheng, Deep learning model with individualized fine-tuning for dynamic and beat-to-beat blood pressure estimation, in: 2021 IEEE 17th International Conference on Wearable and Implantable Body Sensor Networks (BSN), IEEE, 2021, pp. 1–4.
  * [59] P. H. Charlton, J. Mariscal Harana, S. Vennin, Y. Li, P. Chowienczyk, J. Alastruey, Modeling arterial pulse waves in healthy aging: a database for in silico evaluation of hemodynamics and pulse wave indexes, American Journal of Physiology-Heart and Circulatory Physiology 317 (5) (2019) H1062–H1085.
  * [60] I. S. Association, et al., Ieee standard for wearable cuffless blood pressure measuring devices, IEEE Std (2014) 1708–2014.