Machine learning in cardiovascular flows modeling: Predicting pulse wave propagation from non-invasive clinical measurements using physics-informed deep learning

05/13/2019 ∙ by Georgios Kissas, et al. ∙ University of Pennsylvania 0

Advances in computational science offer a principled pipeline for predictive modeling of cardiovascular flows and aspire to provide a valuable tool for monitoring, diagnostics and surgical planning. Such models can be nowadays deployed on large patient-specific topologies of systemic arterial networks and return detailed predictions on flow patterns, wall shear stresses, and pulse wave propagation. However, their success heavily relies on tedious pre-processing and calibration procedures that typically induce a significant computational cost, thus hampering their clinical applicability. In this work we put forth a machine learning framework that enables the seamless synthesis of non-invasive in-vivo measurement techniques and computational flow dynamics models derived from first physical principles. We illustrate this new paradigm by showing how one-dimensional models of pulsatile flow can be used to constrain the output of deep neural networks such that their predictions satisfy the conservation of mass and momentum principles. Once trained on noisy and scattered clinical data of flow and wall displacement, these networks can return physically consistent predictions for velocity, pressure and wall displacement pulse wave propagation, all without the need to employ conventional simulators. A simple post-processing of these outputs can also provide a cheap and effective way for estimating Windkessel model parameters that are required for the calibration of traditional computational models. The effectiveness of the proposed techniques is demonstrated through a series of prototype benchmarks, as well as a realistic clinical case involving in-vivo measurements near the aorta/carotid bifurcation of a healthy human subject.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Recent advances in clinical measurement and computational modeling techniques introduce new capabilities for monitoring the human cardiovascular system from different perspectives such as disease surveys Rose et al. (1982), bio-medical image processing Kak et al. (2002); Haris et al. (2014), computational mathematics Formaggia et al. (2010); Figueroa et al. (2006); Reymond et al. (2009); Grinberg et al. (2011), bio-physics Stefanovska (1999); Ma et al. (2018), etc. These studies reveal the crucial role played by blood flow, arterial wall mechanics and pressure wave propagation, and how their interplay directly characterizes the functionality of the cardiovascular system both in health and disease (e.g., hypertensive disorders O’Rourke (1995)).

Understanding the inner workings of the cardiovascular system has been central to many studies involving clinical, interventional or computational approaches. For instance, Chan et. al. Chan et al. (2007) proposed placing sensors in the human body to achieve real time measuring of arterial velocity and wall displacement to monitor the health of patients. Although the collected in-vivo measurements can be highly accurate, such interventional techniques are sometimes expensive and suffer from limitations that are not easy to address, e.g., difficulties of placing probes in cerebral or uteroplacental arteries Kett-White et al. (2002). These limitations motivate the use of non-invasive measurement techniques such as bio-medical imaging, advances in which currently define the clinical standard of care Edelstein et al. (1980); Antiga et al. (2008); Yushkevich et al. (2006); CIBC (2016). To this end, one of the most commonly used techniques is Doppler ultrasound velocimetry Aaslid et al. (1982) that enables the recovery of blood velocity wave propagation. Using information related to the pulsatility of the flow Mitchell et al. (2011), clinicians are able to infer quantitative information about the pressure in the arterial vessels, e.g., large pulsatility is caused by large pressure Ma et al. (2018). Rather than Doppler ultrasound, more recent studies leverage advances in 4D flow Magnetic Resonance Imaging (MRI) Markl et al. (2012) to recover the full three-dimensional velocity flow field. This method provides a more detailed representation of the flow, and also a larger spatial coverage that can capture quantities like vessel tortuosity in a more accurate manner. Moreover, structural characteristics like the cross-sectional area and vessel diameter can be recovered by segmenting 2D cine images Plein et al. (2001) from MRI measurements. However, critical variables such as the pressure cannot be directly measured by a non-invasive technique Barker et al. (1985), and accurate measurements are only accessible by inserting a catheter equipped with sensors inside the vessel of interest.

These difficulties of directly measuring quantities of interest like the pressure in an in-vivo and non-invasive manner have motivated the use of computer simulations and computational fluid dynamics models to predict them in-silico. Advances in algorithms and computing resources now allow us to perform detailed flow simulations in complex patient-specific arterial topologies using three-dimensional (3D) and/or one-dimensional (1D) formulations of the unsteady Navier-Stokes equations Formaggia et al. (2010); Figueroa et al. (2006); Reymond et al. (2009); Grinberg et al. (2011); Urquiza et al. (2006); Sherwin et al. (2003b). Such tools have been successfully validated against both in-vitro experiments Matthys et al. (2007) as well as in-vivo clinical data Revie et al. (2013), and provide a valuable platform for parametric sensitivity studies Chen et al. (2013). Despite their predictive power, computational models have still not made their way into clinical practice primarily due to their high computational cost and the tedious procedures needed for their practical deployment, e.g., mesh generation, parameter calibration, etc. For instance, such models require the precise subscription of boundary conditions that effectively capture the downstream flow dynamics in small arteries and arterioles via the use of Windkessel or structured tree models Westerhof et al. (2009); Grinberg and Karniadakis (2008); Olufsen (1999); Perdikaris et al. (2015). Inaccurate calibration of the parameters associated with these boundary conditions is often the cause of brittleness in the resulting predictions, thus limiting the translational impact of computational models to the clinical domain.

This goal of this work is to put forth a new paradigm for seamlessly blending non-invasive measurement techniques such 4D flow MRI with computational fluid dynamics models, towards providing a cheap and effective tool for characterizing velocity and pressure pulse wave propagation in human systemic arteries. Leveraging recent advances in physics-informed machine learning Raissi et al. (2019, 2017a); Yang and Perdikaris (2018); Zhu and Zabaras (2018), we put forth an algorithmic framework for producing physically consistent predictions of flow and pressure wave propagation directly from processing noisy and scattered measurements of blood velocity and wall displacement obtained by non-invasive 4D flow MRI. Specifically, we propose to employ deep neural networks to represent the unknown flow variables (blood velocity, wall displacement and pressure) in a given arterial network. We then train these networks to produce outputs that: (i) fit any available clinical data (typically velocity and wall displacement data at a few cross-sections), (ii) satisfy the underlying physical conservation laws as described by a reduced one-dimensional model of pulsatile blood flow Sherwin et al. (2003b); Formaggia et al. (2010), and (iii) ensure conservation and information propagtion across interface points in the network (e.g., bifurcations, junctions, etc.).

In contrast to traditional computational fluid dynamics models, the proposed methodology does not the cumbersome generation of computational meshes, nor it requires the precise prescription of boundary conditions. Instead, it can directly use noisy and scattered measurements obtained by medical imaging at points other than the boundaries to accurately predict the pressure and recover flow information. Furthermore, as flow conditions in general do not vary greatly from one patient to another, pre-trained neural network models can be rapidly adapted to a new patient case. After the model is trained, all quantities of interest (blood velocity, wall displacement and pressure) can be predicted for every temporal or spatial point in the network. This directly enables the calculation of Windkessel parameters as a simple post-processing step, leading to a simple scheme for calibrating more detailed and expensive 3D models of blood flow.

The rest of the paper is organized as follows. In section (2), we provide an outline of the one-dimensional blood flow model employed in this work, as well as the deep learning techniques used for obtaining physically consistent predictions of pulse wave propagation directly from clinical measurements. In section (3), we demonstrate the effectiveness of the proposed methods in a series of synthetic systematic studies, as well as a realistic clinical case involving in-vivo measurements near the aorta/carotid bifurcation of a healthy human subject. Conclusions and further discussions are given in section (4). All code and data accompanying this manuscript will be made publicly available at https://github.com/PredictiveIntelligenceLab/1DBloodFlowPINNs.

2 Methods

2.1 A simplified model of pulsatile flow in arterial networks

Figure 1: A schematic representation of a simple bifurcating arterial topology and its one-dimensional center-lines used for the reduced model. Throughout this work, three dimensional geometries are recovered from MRI data and the corresponding center-lines are extracted by using the VMTK library Antiga et al. (2008). Blood velocity data are typically obtained using Doppler ultrasound or 4D flow MRI, while the area data are parsed from 2D Cine images recovered by 4D flow MRI.

Pulse wave propagation in arterial networks can be effectively modeled using one-dimensional (1D) reduced order models Sherwin et al. (2003b); Reymond et al. (2009, 2011). In order to achieve the order reduction, a series of assumptions need to be made Sherwin et al. (2003a). First, we assume that the local curvature is small enough such that the geometry can be described using a Cartesian coordinate , as shown in figure 1. Moreover, the fluid is incompressible and Newtonian, since we are considering geometries consisting of large arteries, so the density and dynamic viscosity are constant. Lastly, the structural properties of the artery are preserved at a cross-section. Following Sherwin et al. (2003a); Formaggia et al. (2010) we consider a reduced form of the incompressible Navier-Stokes equations. Conservation of mass and momentum can be expressed as a hyperbolic conservation law Sherwin et al. (2003a) that describes the evolution of blood velocity and cross-sectional area Sherwin et al. (2003b); Formaggia et al. (2010). In order to close the system, a third equation accounting for the relation between the pressure and the area is used, which is derived by assuming a thin wall tube and using Laplace’s law Sherwin et al. (2003a). This reduced order model provides an accurate representation of the underlying transport processes and its effectiveness in correctly capturing pulse wave propagation phenomena has been validated against both in-vitro and in-vivo data Matthys et al. (2007); Reymond et al. (2009, 2011).

The system derived by the above analysis takes the form Sherwin et al. (2003b); Formaggia et al. (2010)

(1)

where , and represent the cross-sectional area, the velocity and pressure, respectively. In addition, and represent the spatial and temporal coordinate in each vessel, respectively. Also, and denotes the vessel’s sectional area at equilibrium, the wall thickness, the Young’s modulus, the external pressure, the Poisson ratio and is the density of the blood. Values for , and are typically mined from the literature, while and are typically set to the diastolic pressure and wall displacement, respectively. For the examples presented in section (3) the parameter is computed based on the empirical relation in Olufsen (1999). Without loss of generality, we will not take into consideration wall viscous effects.

2.2 Physics-informed neural networks

Physics informed neural networks (PINNs) is a new paradigm that leverages recent advances in deep learning to infer solutions, parameters and/or constitutive laws invlolving partial differential equations (PDEs)

Raissi et al. (2017b, c, 2019); Tartakovsky et al. (2018). In this framework, the solution of partial differential equations is parametrized by a neural network that is trained to match the measurements of the system, while being constrained to approximately satisfy the underlying physical laws. In our particular case, we define one neural network to represent the solution of the equation (1) for vessel in our arterial network. Here denotes the parameters of the neural network for vessel . The network aims to approximate the following mapping:

In correspondence to equation (1), the following residuals can be defined:

(2)

These residuals can be then used as constrains during the training of the neural networks in order to encourage them to produce physically consistent predictions. In addition to minimizing the residuals, the neural networks are also constrained to fit any available non-invasive clinical measurements that may be noisy and scattered in space and time. In order to illustrate our model, we refer figure (2). For example, in each vessel , we have scattered measurements of (corresponding to the black ’s in figure (2)), namely , , and (corresponding to the red ’s in figure (2)), namely , , along with a larger number of collocation points , , that aim to satisfy the constraints at a finite set of collocation nodes. In the examples presented in this paper, we want to emphasize that we will only use measurements for and , but not for

, as these are the quantities that often cannot be reliably measured non-invasively in the clinic. More details on how to construct the respective terms of the loss function for training the neural networks

will be given in the next section.

Figure 2: Graph of a prototype arterial network where 1 vessel splits into 2 vessels at the bifurcation point: (a) Structure of the simple arterial network consisting of 3 vessels with 1 bifurcation point. (b) Domain of vessel . (c) Domain of vessel . (d) Domain of vessel . Blue points represent the collocation points. Green stars represent the interface points at the bifurcation. Red crosses indicate the locations of the velocity measurements. Black crosses show the locations for the area measurements.

2.3 Loss function: Measurements

This part of the loss function corresponds to fitting the clinical measurements obtained for some of the vessels in the network. In figure (2) these are denoted as black and red ’s, the color depending on the quantity we have measurement for. In this particular case, this term encourages the output of the neural network and to match the measurements of area and velocity obtained by a clinical procedure (e.g., segmenting 2D cine images and Doppler ultrasound Aaslid et al. (1982)). This part of the loss function has the form (take vessel as example):

(3)

where , represent the number of measurements for and in vessel respectively. Also, denote total number of vessels in which we have measurements. In the above equation and represent the outputs given the parameters of the neural network for vessel . Minimizing this term will encourage the neural networks to fit the available measurements.

2.4 Loss function: Collocation points

This part of the loss function corresponds to the collocation points, see the blue dots in figure (2). These are points that we randomly choose inside the arterial domains using a latin-hypercube sampling strategy Stein (1987). Over these collocation points, we impose the physical constraints by encouraging the right hand side of equation (2) to be equal to zero. The partial derivatives in the residual expression (equation (2)) can be computed using automatic differentiation Paszke et al. (2017); Abadi et al. (2016). This objective is imposed to encourage the neural networks to find a particular set of parameters that make their predictions consistent with the underlying differential equations, which translates to having the minimum residual. By satisfying this condition together with matching the measurements at particular points in section (2.3), we can obtain a model that is capable of inferring the solution at any spatial coordinate of the domain and any time. This collocation loss function takes the following form (take vessel as example):

(4)

where represent the number of collocation points in vessel , meaning many blue dots we have in figure (2). Also, denote the total number of vessels in our arterial network. The terms , and represent the residual of area, velocity and pressure defined in equation (2), respectively.

2.5 Loss function: Interfaces

One-dimensional models can be extended to treat splitting and merging arteries by imposing proper boundary conditions Sherwin et al. (2003a). In this case, one needs to provide boundary conditions for each artery at the interface points to ensure conservation. In conventional numerical methods (e.g., Discontinuous Galerkin method Cockburn et al. (2012); Sherwin et al. (2003b)), the velocity and the area can be discontinuous at the interface points (e.g., bifurcations, junctions), so in order to find the values, a Riemann problem has to be solved. This is typically done by employing the characteristic variables of the hyperbolic system of equation (1), accounting for a decoupled system of scalar equations, so that the travelling waves can reach the splitting points Sherwin et al. (2003b). The proposed method can work without using information on the characteristics, instead just imposing the momentum and mass conservation equations suffices. To illustrate our workflow, let us consider the case where the artery splits to and , see figure (2). For each vessel , the area, velocity, pressure and density are denoted as , the spatial and temporal variables are denoted as . In reference to figure (2), this corresponds to the green stars and the index on how many of these we use to impose the continuity. The boundary conditions at the bifurcation points are derived by the conservation of momentum and mass Sherwin et al. (2003a) for the parent and daughter vessels:

(5)

In the case of more than one bifurcations, these conditions are imposed to every splitting point using the appropriate notation, depending on the numbering of the domains. In the case of junctions, instead of bifurcations, or mixed change of geometry, system (equation 5) takes the form of the appropriate condition following the conservation laws. In general, the interface loss has the form (take bifurcation point as example):

(6)
(7)
(8)

where the indices , , in denote the father and daughter vessels, respectively, at each bifurcation. Also, denotes the total number of bifurcation points in our arterial network. represent the number of collocation points on the interface boundaries, i.e. how many green stars we have in figure (2). So, in the above equation if we choose for example that means we calculate the pressure at the interface point with index using the network corresponding to the father vessel, in the notion of figure (2) this would correspond to domain . For the interface, we feed the neural network the inputs , where is the coordinate of the interface point. Minimizing the interface loss encourages the neural network to satisfy the conservation laws at the interface points.

The interface loss function ensures the information of the measurements in one domain can be propagated throughout the neighboring domains. In the original formulation of Sherwin et. al. Sherwin et al. (2003a) a Riemann problem was solved on the interfaces due to the hyperbolic nature of the wave propagation equations. Here we have experimented with including a constraint ensuring the propagation of characteristic variables at interfaces, but we did not observe any substantial improvement in terms of accuracy and computational efficiency, thus we choose to not include these additional constraints in the neural networks’ loss function.

2.6 Loss function: Accounting for all constraints

In the previous sections, we provided details on how to construct the individual terms of the loss function corresponding to the measurements, physical constraints and continuity, that the model should satisfy. In this section, we will show how to combine these individual terms to create the joint form of the loss function, together with some examples to clarify the procedure. The loss function of our physics informed neural networks (PINNs) takes the form:

(9)

according to the definition in previous section, represent the total number of domains, denotes the indices of domains where we have data and correspond to the number of interfaces in the arterial networks. The index denotes the identity of the domain under consideration and the index is helping as keep track of the interface points. Index does not refer to specific neural network, because , and on the bifurcation points are calculated through the father and daughter vessels by their respective network. This complete form of the loss function that encourages the model to learn the underlining physics of the problem by using measurements inside the domain, not necessarily at the boundaries e.g., figure (2), and by respecting the imposed physical constraints. We will present how the loss function is constructed for some simple cases in order for the method to be more clear to the reader.

Figure 3: Graph of a complex arterial network which involves 7 vessels and 3 bifurcation points: The red points represent the boundaries of the arterial networks. The greed points represent the bifurcation points where we impose the continuity boundary conditions. The number on each vessels are the corresponding label for these vessels. The vertical lines denote the middle points of the vessels and also the spatial locations for which we have training data, e.g., for the velocity (red) and the wall displacement (black).

For a single artery, the loss function for our method is the following:

(10)

In this particular geometry there are no bifurcations, so , thus the term corresponding to splitting, , is neglected.

For the case that we have three vessels and one bifurcation, see figure (2), we specify the loss function as:

(11)

where the first summation corresponds to the domains where we have measurements for the velocity and cross-sectional area. The second summation encourages the model to satisfy the continuity equation on the bifurcation point. The third summation represents the physics-informed constraints in all three vessels. The second term is written as a summation, for consistency, but in this case, , so there is only one component of .

For a network with 7 vessels and 3 bifurcations, see figure (3), we specify the loss function as:

(12)

Here, the first summation corresponds to the seven vessels where we have measurements of the velocity and area. The second summation correspond to the three bifurcation points. The third summation is representing the physics-informed constraints on seven different arteries. The loss is minimized by encouraging the model to reconstruct the measurements provided at certain points where the measurements are obtained and to respect the constraints imposed by the physics of the problem, e.g., the differential and continuity equations. The loss function in equation (12

) is minimized using stochastic gradient descent algorithm

Ruder (2016) or its modern variants Kingma and Ba (2014); Duchi et al. (2011) that are widely used in current deep learning research.

2.7 Non-dimensionalization and normalization

In equation (1), the order of magnitude of the different physical quantities, velocity, cross-sectional area and pressure, have a significant relative difference, e.g., , and , which casts great difficulty on the training of the neural network Glorot and Bengio (2010). A direct application of the methods proposed in Raissi et al. (2017b, c, 2019) cannot handle this issue. For overcoming this problem, we employ a non-dimensionalization and normalization technique with the purpose of scaling the input and the output of the neural networks in a proper scale (e.g.,

) and normalizing the spatial and temporal coordinates to have zero mean and unit variance for training the neural networks more efficiently. For the purpose of non-dimensionalization we introduce the characteristic variables, which are commonly used in multi-scale physics

Famiglietti and Wood (1994) in order to simplify the equations. For this problem we need a characteristic length and a characteristic velocity .

We will choose the characteristic length to be the square root of the mean of the equilibrium cross-sectional area of the network vessels. In order to choose the characteristic velocity we make use of the physiological condition that the wave speed in a vessel has to be one order of magnitude larger than the length Sherwin et al. (2003b), given that, in the normalized length case . Thus:

where . At this point we define the quantities:

(13)

where , and .

Now, the input of the neural network has support . It is shown that normalizing the input to have zero mean and unit variance makes the training of the neural network more efficient as it prevents gradient saturation and provides stable updates Glorot and Bengio (2010). In this step, we will apply this technique to and , in the vessel , and get:

(14)

where , the mean value and ,

the standard deviation of the spatial and temporal coordinates for the vessel

, respectively, and , the scaled inputs. Using the variables stated above, we derive the updated system of equations (take vessel as example):

(15)

In this non-dimensional and normalized form, all the variables and inputs are scaled to order . This is what effectively enables the training of our physics-informed neural networks in this complex setting. Likewise, at the predicting stage, we first scale the inputs by the characteristic variables as equation (13), i.e. by , then normalize them following equation (14). Finally, we revert the predicted quantities back to their original form by multiplying with the characteristic variables i.e. .

In order to be consistent with the derivation above, we have to follow the same procedure for every condition we impose to the model. In this notion we derive the non-dimensional continuity conditions, by inserting the non-dimensionalizing quantities into the conservation laws. By doing so, we get:

(16)
(17)
(18)

The above system of equations must be satisfied at any interface, which is achieved by the second term in the loss function in equation (9) that forces the model to respect the conservation laws. Finally, the non-dimensional equations (15), (16), (17) and (18) define the optimization objectives that are used in equation (9) for minimizing the residual at the collocation, bifurcation and measurement points respectively. As mentioned above, we have to multiply the predictions of the network by the scaling parameters , and , when doing inference or else we get a scaled version of the solutions.

2.8 Windkessel model parameter identification as post processing

A crucial aspect in calibrating computational models of blood flow is related to the proper choice of outflow boundary conditions. Due to the excessive computational cost of imaging and performing pulsatile flow simulation across the entire human arterial tree, typically only truncated topologies are considered. This introduces the need of modeling the neglected downstream dynamics in smaller arteries and arterioles; a task that is typically accomplished using lumped parameter Windkessel models Westerhof et al. (2009). Windkessel models provide us with a way to account for arterial blood pressure in terms of the compliance of the large elastic arteries and the resistance of the smaller arteries, by making an analogy with an electric circuit. In this setup, the compliance plays the role of a capacitor and the resistance of a resistor and the whole circuit helps to correctly model the influence of the small downstream vessels to the blood circulation. In our case we consider the popular three-element Windkessel model with governing equation (take vessel as example):

(19)

where represent the total number of outlets where we want to identify the parameters. In the equation above, is the flow rate at the outlet, Pa denote the downstream pressure that is constant for all vessels. Also, represents the characteristic impedance Alastruey et al. (2008) that can be easily computed using , the total arterial compliance and the systemic vascular resistance. The characteristic impedance is chosen in a way that allows the incoming wave to reach and without being reflected Alastruey et al. (2008). The main challenge here relates to choosing the compliance and resistance parameters and such that physiologically correct results are sought.

Some of the work done in estimating Windkessel model parameters is based on modeling assumptions. Grinberg et. al. Grinberg and Karniadakis (2008) proposed a method of imposing a simple resistance boundary condition using the ratio of flow rates of terminal vessels. They also showed that this method is robust and stable and scales to large arterial systems. The accuracy can be guaranteed by choosing appropriate values for the resistance at each outlet. Ismail et. al. Ismail et al. (2013) proposed a method for tuning the parameters of a three-element Windkessel parameters for a coupled 3D/0D system by implementing an inverse analysis based on adjoint methods, using the patient-specific pressure information as a target. Spilker et. al. Spilker and Taylor (2010) proposed a multidomain method that gets coupled with a reduced order, a three-element Windkessel, model in which by imposing characteristic information of input pressure, one can tune hemodynamic simulations. The pressure characteristics, such as minimum or maximum pressure, are imposed in the input of the domain. Also Spilker et. al. Spilker et al. (2007) proposed a procedure for imposing and tuning impedance type boundary conditions by adjusting the length of morphometry produced trees to match the outflow. Yih-Choung Yu et. al. Yu et al. (1998)

presented an online estimator based on extended Kalman filter that uses sets of measurements to estimate state variables and at the same time recover the model parameters. Pant

et. al. Pant et al. (2014) presented an algorithm that iterates between a three dimensional and a zero-dimensional model with known pressure, at specific points, to recover the values of three-elemnent Windkessel model using Unscented Kalman filters. Sciavazzi et. al. Schiavazzi et al. (2017) proposed a method for automatically estimating Norwood circulation model parameters and their uncertainty related to the clinical measurements. In this framework an analysis over the confidence of the estimations can be performed using a differential evolution adaptive Metropolis estimation and patient specific measurements on pressure and flow.

All the aforementioned works require the precise prescription of inflow and boundary conditions and rely on the repeated evaluation of conventional solvers for different combinations of and until a parametric setting that yields a reasonable match to clinical data is identified. For networks with more than a handful of outlets this defines a tedious and prohibitively costly calibration procedure. In contrast, our proposed methodology can yield a cheap and effective procedure for calibrating Windkessel model parameters without the need of employing conventional simulators. Specifically, the training of the proposed physics-informed neural networks (as described in section (2)) does not require the a-priori specification of Windkessel-type outflow boundary conditions. Moreover, once the neural networks are trained on scattered velocity and wall displacement data, they can yield physically consistent predictions for the propagating pressure wave. This information can be then used for identifying the resistance and compliance parameters and at each outlet by solving the differential equation (19) as a post-processing step.

To this end, here we propose a simple method for identifying the parameters and in equation (19) by adaptive grid search given the predicted outflow data and produced by the physics-informed neural networks representing the flow solution at each outlet. Once the model is trained we can predict and at every locations in the domain and at any times. Thus we have access to the data as a time series and via the neural network output, where is the length of the time series. We select a periodic part of the data and fit using Fourier series consisting of fifty modes and then compute the time derivative of with respect to time . Thus, we are setting up a post processing step that leverages the solution for and its derivative at each temporal discretization and use the ODE solver to compute using some combination of parameters . In this case, we compare the difference between and the pressure predicted by the model. Adaptive mesh refinement of the parameter space can be used to discover the optimal solution for by iteratively exploring new meshes after computing the loss function. We define the loss function as a function of at points randomly selected within the interval :

(20)

In this case, we compute at the parameter grid and find the set of values for which this loss gets minimized. Then we explore an new interval , and refine the mesh in this vicinity. This is done for a number of consecutive times. In our case, we consider the number of consecutive times as 5, in order to reach the parameter set that provides us with the smallest value of loss in equation (20). The computational time required for discovering the Windkessel parameters is about 10 minutes per outlet using a conventional computer. This interval could vary depending on the number of iterations and the number of solution points provided. This is just a fraction of the time required to train the network and a computationally cheap procedure compared to other methods, e.g., adjoint optimization Ismail et al. (2013).

3 Results

The proposed method will be tested in three different cases. In the first case, we consider a prototype arterial network resembling a carotid bifurcation that consists of 3 vessels with 1 bifurcation as shown in figure (2). In the second case, we consider a complex arterial network resembling an idealized pelvic geometry that consists of 7 vessels and 3 bifurcations as shown in figure (3). For both of these two cases the data used to train the model are synthetic data generated using a conventional Discontinuous Galerkin simulator Sherwin et al. (2003b)

. More specifically, the data set is created by considering the results of the Discontinuous Galerkin method for velocity and wall displacement at a steady state cardiac cycle. For the first case, we also provide a series of comprehensive systematic studies in order to quantify the accuracy and robustness of the proposed methods. Finally, we also present results for a real clinical case utilizing 4D flow MRI data collected near the aorta/carotid bifurcation of a healthy human subject. The proposed algorithms are implemented in Tensorflow v1.10

Abadi et al. (2016), and computations were performed in single precision arithmetic on a single NVIDIA Tesla P100 GPU card. The neural network models for the first and the third case was trained for, approximately, 1 hour and 20 minutes and for the second case for approximately 7 hours. All code and data accompanying this manuscript will be made publicly available at https://github.com/PredictiveIntelligenceLab/1DBloodFlowPINNs.

3.1 Flow through a prototype Y-shaped bifurcation

3.1.1 Arterial network topology and observed measurements

Here we consider a simple arterial network resembling a prototype carotid bifurcation. The network consists of 3 vessels with 1 bifurcation where each vessel has lengths mm, mm, mm, respectively, with one bifurcation point at mm. Somewhere within these vessels we assume we have access to time-series data for the blood velocity and the wall displacement, as is shown in figure (2). Specifically, we assume measurements at the inlet and outlets of the arterial network for which represents approximately four cardiac cycles with a period of . We should emphasize that we do not assume any measurements for the pressure at any location in the arterial network.

3.1.2 Generation of synthetic training and validation data

For this example we created a set of synthetic data for the cross sectional area and velocity using an in-house Discontinuous Galerkin (DG) solver Sherwin et al. (2003b); Perdikaris and Karniadakis (2014); Perdikaris et al. (2015). Out of the resulting waveforms we choose 4 cycles from the steady state solution as training data for the neural network model. For the DG simulation, we considered nominal values for the blood density Kg/m and the blood viscosity mPas, respectively Fossan et al. (2018). Moreover we provide the DG solver with the Windkessel and arterial wall parameters shown in table (1).

Arterial Length Peripheral Peripheral Equilibrium cross-
segment (cm) Resistance compliance sectional area
(10 Pa s m) (10 m Pa) (Pa/m) (m)
1 17.03 - - 6.97e+07 1.36e-05
2 0.7 0.5251 0.3428 5.42e+08 1.81e-06
3 0.67 0.2702 0.6661 6.96e+07 1.36e-05
Table 1: Generation of synthetic training and validation data for the prototype Y-shaped bifurcation: Physiological data used as simulation parameters for the Discontinuous Galerkin method.

3.1.3 Neural network approximation set-up

We parametrize the solution of this problem using three neural networks, one for each vessel. Each of these networks consists of 7 hidden layers with 100 neurons per layer, followed by hyperbolic tangent activation function. The neural networks are intiialized using Xavier initialization

(Glorot and Bengio, 2010). This architecture has enough capacity to efficiently capture fine features in the propagating waveforms. A more detailed discussion on the neural network structure can be found in section (3.1.5) where we present a systematic study that justifies the aforementioned choices in terms of balancing the trade-off of computational cost and predictive accuracy. The neural networks take the scaled inputs and predict the non-dimensional outputs , as discussed in section (2.7). We are providing each neural network with velocity and cross-sectional area time series for one point at the inlet and the two outlets of the arterial network. Using this data we infer the solution at , randomly chosen (using Latin hypercude sampling Stein (1987)) collocation points, inside each domain, and interface points. Moreover we randomly choose batches of points as input to the Adam optimizer Kingma and Ba (2014) and learning rate for the first iterations. Consequently, we set for the following for further minimizing the error and fine-tuning the neural network parameters.

3.1.4 Comparison of model predictions and reference solution

We first demonstrate the importance of the non-dimensionalization and normalization of the governing equations, as described in section (2.7). In order to compute the unknown quantities we assume that we have access to the values of and for all times at one spatial point in each vessel, figure (2). We show in figure (4) that the original methodology introduced in Raissi et al. (2017b, c, 2019) fails to overcome the difficulty of handling the different scales and thus the numerical prediction collapses. As we mentioned in section (2.7), this phenomenon occurs due to the significant difference between the orders of magnitude of the model variables, i.e. , , and , which casts great difficulty on the training of the neural network. We, also, show that the proposed normalization and non-dimensionalization tools significantly improve the model training capabilities and, thus, provide accurate predictions.

Figure 4: Flow through a prototype Y-shaped bifurcation: (a) Comparison of predicted velocity wave between Discontinuous Galerkin (red), physics-informed neural networks with non-dimensionalization (blue) and physics informed neural networks without non-dimensionalization (black) at the middle point of vessel . (b) Comparison of predicted pressure wave between Discontinuous Galerkin (red), physics-informed neural networks with non-dimensionalization (blue) and physics-informed neural networks without non-dimensionalization (black) at the middle point of vessel .

The results produced by our model for the middle points of each domain and randomly selected values of are presented in figure (5). Evidently, there is a good match between the solution acquired from Discontinuous Galerkin solver Cockburn et al. (2012) and the one predicted by the proposed physic-informed neural networks model. The prediction error is further quantified in the systematic studies presented in section (3.1.5).

Figure 5: Flow through a prototype Y-shaped bifurcation: Comparison of the physics-informed neural network model predictions and the discontinuous Galerkin reference solution. The numbering of the vessels starts from left to right, with the far left figure corresponding to the father (vessel ) and the middle and right one to the daughter (vessels and , respectively), as shown in figure (2). The top row represents the predicted velocity and the bottom row the predicted pressure for each vessel.

To highlight the role of the continuity equation described in section (2.5), we present a comparison between the left hand side and the right hand side of the continuity equations (16) and (17),(18) in figure (6). We can infer from the figure that the mass in vessel is equal to the total mass in vessel and vessel . Moreover, the momentum at the father and daughter vessels at the bifurcation point is equal to each other. Thus, the conservation laws are well preserved in our model.

Figure 6: Flow through a prototype Y-shaped bifurcation: Comparison of conserved quantities at the bifurcation point after training the proposed physics-informed neural networks. (a) Conservation of mass (see equation 16). (b) Conservation of momentum (see equations 17 and 18) for the bifurcation point.

3.1.5 Systematic studies

We present a series of comprehensive systematic studies to test the sensitivity of the proposed methods and quantify their robustness with respect to different hyper-parameter settings. All numerical experiments are performed for the canonical Y-shaped bifurcation topology presented in the previous section. To test the sensitivity of the proposed methods with respect to the initialization of the neural networks, we consider a data set comprising of , , , training, collocation, interface points and batch size, respectively, and fix the architecture for our neural networks to consist of seven hidden layers containing hundred neurons each followed by a hyperbolic tangent activation function. By doing so, we create an ensemble of fifty cases all starting from a Xavier initialization (Glorot and Bengio, 2010) for all network weights (with a randomized seed), and a zero initialization for all bias parameters. In table (2) we report the relative error between the predicted mean solution and the known Discontinuous Galerkin solution Cockburn et al. (2012) for all 50 randomized trials at point in vessel for randomly selected temporal values. Evidently, our results are robust with respect to the the neural network initialization as in all cases the stochastic gradient descent Kingma and Ba (2014) training procedure converged to solutions that are close to each other. We summarize the result by reporting the mean and the standard deviation of the relative error as

Relative error
3.42e-02 4.99e-02 5.65e-02 3.92e-02 1.07e-01 5.52e-02 5.34e-02 3.54e-02 6.41e-02 4.57e-02
4.62e-02 4.04e-02 2.94e-02 3.78e-02 3.48e-02 1.32e-01 6.01e-02 3.00e-02 4.78e-02 4.63e-02
3.68e-02 4.48e-02 4.37e-02 5.81e-02 1.23e-01 5.80e-02 4.06e-02 2.33e-02 3.02e-02 3.96e-02
5.70e-02 3.61e-02 3.85e-02 3.83e-02 4.82e-02 4.85e-02 5.46e-02 6.01e-02 3.84e-02 7.43e-02
5.82e-02 5.09e-02 4.13e-02 6.66e-02 1.31e-01 7.31e-02 6.09e-02 5.59e-02 3.50e-02 3.47e-02
Table 2: Systematic study on the neural network initialization: Relative prediction error for different neural network initialization using different randomized seeds. The prediction errors are obtained by comparing the predicted value of pressure with the reference one in parent vessel at point for randomly selected temporal values.

Subsequently, we test the sensitivity of our model with respect to the architecture of the neural networks. In this case, we fix the number of noise-free training data as , and for each vessel. In all cases, we used a hyperbolic tangent activation function. We initialize our neural network weights using Xavier initialization Glorot and Bengio (2010). In table (3) we report the relative prediction error for different fully connected neural network architectures (i.e. different number of layers and number of nodes in each layer). The general trend suggests that as the neural network capacity is increased, we obtain more accurate predictions, which is indicating that our physics-informed constraint on the PDEs residual can effectively regularize the training process and safe-guard against over-fitting. For the rest of the paper we will utilize an architecture defined by seven layers and a hundred neuron per layer for its good performance in balancing the accuracy and computational intensity.

20 50 100 200
1 6.14e-01 8.41e-01 7.04e-01 1.07e+00
3 4.09e-01 3.08e-01 3.68e-01 2.39e-01
5 5.18e-01 2.35e-01 1.03e-01 2.49e-01
7 4.42e-01 3.88e-02 3.73e-02 3.56e-02
Table 3: Systematic study on the neural network architecture: Relative prediction error for different fully connected neural network architectures with different number of hidden layers , and neurons in each layer. The prediction errors are obtained by comparing the predicted value of pressure with the reference one in parent vessel at point for randomly selected temporal values.

3.1.6 Windkessel parameter identification by adaptive grid search

Figure 7: RCR identification by adaptive grid search for a Y-shaped bifurcation: (a) A sketch that illustrates the adaptive grid search. (b) The surface illustrating the values of the loss with respect to parameters and . The red star denotes the point in the parameter space where the loss function achieves the minimum value.

As described in section (2.8), our method provides an effective mechanism for c//;/alibrating Windkessel model parameters as a simple post-processing step. Specifically, using the predicted flow and pressure at the outflow points, we start by defining a coarse mesh in the space of Windkessel model parameters, and use an ODE solver with inputs , and the inferred flow to calculate the pressure wave form. After we calculate the pressure waveform for every grid point, we find the one that produces the smallest relative error when compared with the predicted one and refine the mesh near its vicinity. By repeating this procedure for a number of times, in this case five, we acquire a set of and that produce the minimum relative error when compared with the target one. In figure (7) we present the relative error as a function of the parameters and for the case of vessel as a visual example of how the method of parameters identification works. For details about the predicted , and the corresponding relative loss, equation (20), we refer to table (4).

vesselParam P waveforms
7.58e+08 4.89e-11 3.36e-02
1.84e+09 7.36e-11 3.67e-02
Table 4: Windkessel parameter identification for a Y-shaped bifurcation: Resistance and compliance discovered by adaptive grid search for the three vessels, one bifurcation case and the relative error between the pressure waveforms predicted by the physics-informed neural networks and the one reconstructed using equation 19 with the identified and parameters.

3.2 Flow through an idealized pelvic arterial network

3.2.1 Network topology and observed measurements

Our goal here is to demonstrate the effectiveness of the proposed methodology for a slightly more complicated problem, where the available training data lies inside the arterial domain, i.e., no available data is assumed at the inflow or the outflow boundaries. This illustrates a favorable property of our approach versus traditional simulation tools that necessitate the precise prescription of boundary conditions, especially for cases in which measurements are not accessible on the boundaries. To this end, we consider a network consisting of seven arteries with three bifurcations, where every single artery splits into two. This topology resembles an idealized female pelvic geometry, starting from the descending aorta, splitting down to the uterine arteries. We set the vessel lengths to be , , , , , and , as shown in figure (3). Each vessel is considered to have its own local coordinate system defined in the interval , where is the length of the respective vessel. The arterial topology includes three bifurcations at the end points of domains , and , as depicted in figure (3). For this geometry we assume we obtained measurements for the cross-sectional area and velocity over time at the middle points of domains , , , and and we want to predict the pressure at the boundary points. No information is provided for domains and , and no information about the pressure in any vessel. In order to show that in this case the model can also propagate the information via the boundaries, we aim to infer the pressure at two points in domains and namely and .

3.2.2 Generation of synthetic training and validation data

In this example we generated a set of synthetic data for cross-sectional area and velocity using an in-house Discontinous Galerkin solver Sherwin et al. (2003b); Perdikaris and Karniadakis (2014); Perdikaris et al. (2015). Out of the resulting waveforms, we chose 3 steady state cycles as training data. For the DG simulation, we choose the blood density to be equal to Kg/m and the viscosity mPas Fossan et al. (2018). Moreover we provide the DG solver with the Windkessel and structural parameters shown in table (5).

Arterial Length Peripheral Peripheral Equilibrium cross-
segment (cm) resistance compliance sectional area
(10 Pa s m) (10 m Pa) (Pa/m) (m)
1 1.0682 - - 2.65e+07 2.14e-05
2 6.66638 - - 2.60e+07 2.21e-05
3 6.99352 - - 2.63e+07 2.17e-05
4 14.7735 0.3133 16.62 2.82e+07 1.97e-05
5 14.9503 0.1654 31.49 2.71e+07 2.08e-05
6 13.6421 0.1682 30.96 2.67e+07 2.12e-05
7 13.4384 0.2086 2.092 2.87e+07 1.92e-05
Table 5: Flow through an idealized pelvic arterial network: Physiological data used as simulation parameters for the Discontinuous Galerkin method.

3.2.3 Neural network approximation set-up

We parametrize the solution of each artery by a neural network, seven in total, consisting of seven hidden layer and one hundred neurons followed by a hyperbolic tangent activation function, each. The neural network takes the scaled inputs and predicts the non-dimensional outputs for each domain. For domains , we provide only the initial conditions, and no other information on and , and train the model to be able to predict these values inside these domain by employing information propagated by the conservation laws on the boundaries. We are randomly, using Latin-hypercube sampling Stein (1987), select collocation points in each domain and have measurements of and , sampled at the middle points of domains , , , and . For domains and we provide as inputs, where the subscript denotes the coordinates and the time of the initial conditions. In this case is a constant array having zero values and a constant array comprised of values of equispaced points inside the domain. We choose randomly selected points at each training step and the same number of interface points for imposing the continuity conditions. The model is trained employing the Adam optimizer Kingma and Ba (2014) with learning rate for the first iterations and for the consequent . We observe from the results presented in figure (8) that there exists some discrepancy between the predicted values for vessels and for which no data is provided.

3.2.4 Numerical comparison and Windkessel parameter identification

In this case, we aim to calibrate a three-element Windkessel model parameters by employing our algorithm to predict the flow at the outflows of the arterial network (red points in figure (3)) and then utilize these values, as explained in section (2.8), to solve equation (19). We start by creating a coarse mesh for and and solve the ODE on the grid points using the predicted flow. As long as, we calculate the pressure for all grid points, we choose the one for which the relative error between the ODE solution and the model prediction exhibits the smallest value. Then, we generate an identical mesh around its vicinity. As soon as we calculate the waveforms and their respective relative error for the new grid points, we refine the mesh again. We repeat this process for a number of iterations, five in our case, to acquire the optimal values of resistance and compliance. For details on the predicted , and the corresponding relative loss we refer the reader to table (6). In figure (8) the numerical results of our method are presented, in comparison to the solution produced by the discontinuous Galerkin method and the discontinuous Galerkin method using the discovered Windkessel parameters.

Figure 8: Numerical comparison for an idealized pelvic arterial network: On the left the comparison between the predicted velocity and on the right the predicted pressure of the three methods is presented. The comparison is between our model (blue), the reference discontinuous Galerkin solution (red) and the discontinuous Galerkin solution using the discovered Windkessel parameters (black). On the right, the same case is presented. The vessels are numbered form top to bottom, starting with 1 at the top to at the bottom.
vesselParam P waveforms
1.36e+09 3.78e-09 7.92e-03
1.64e+09 2.22e-10 3.29e-02
1.22e+09 2.96e-09 8.26e-03
2.69e+09 1.84e-09 2.56e-02
Table 6: Windkessel parameter identification for an idealized pelvic arterial network: Resistance and compliance discovered by the adaptive grid search for the seven vessels, three bifurcations case and the

relative error between the pressure waveforms of the ordinary differential equation solver using the identified parameters and the model prediction.

3.3 Flow through the aorta/carotid bifurcation of a healthy human subject

For the last example we will test the effectiveness of the proposed method on a realistic clinical case involving measurements acquired in the aorta/carotid bifurcation of a healthy volunteer.

3.3.1 Geometry and problem setup

In this case we consider an arterial geometry consisting of 4 vessels containing 1 bifurcation. An MRI image of this geometry is illustrated in figure (9). By utilizing medical image processing and blood flow velocity measuring techniques (see section (3.3.2)), we acquire measurements of area and velocity at the Aorta , Aorta , Aorta , Aorta and carotid points. Our goal is to leverage the data at the carotid, Aorta and Aorta points to predict the area, velocity and pressure anywhere inside the arterial network. More precisely, we will test the predictions of our model against the measurements at Aorta . Furthermore, we would like to utilize the predicted pressure at the carotid point and Aorta to discover and parameters, as discussed in section (2.8). For vessel we provide the network with just the initial conditions of velocity and area. We adopt a 4 vessels geometry instead of 3, as the equilibrium cross-sectional area is not constant at the vessel ends. Thus, by considering 4 vessels and to be a linear function, we increase the accuracy of the model. For this case, we choose the neural network architecture and parameters to be the same as in section (3.1).

Figure 9: Flow through the aorta/carotid bifurcation of a healthy human subject: Positions of acquired 4D flow MRI measurements in the aorta/carotid bifurcation of a healthy volunteer. Path length measurements are presented in reference to “Aorta1”. Measurements collected at locations [“Aorta1”,“Aorta4”,“Left common carotid artery”] were used to train the physics-informed neural network model, while measurements at locations [“Aorta2”,“Aorta3”] were used to validate its predictions.

3.3.2 Clinical data acquisition

The measurements used in this section were based on acquired MR images from a healthy female volunteer (age=27 years, weight=51 kg, height=160 cm) using a 1.5T Avanto scanner (Siemens Healthcare, Erlangen, Germany). The patient was lying head-first, supine on the table. For the 3-lead electrocardiogram (ECG) monitoring, a neck coil, two body array coils on the chest and abdomen, and a spine array coil were used.The MRI protocol consisted of a Balanced Steady-State Free Precession (bSSFP) localizer of the neck through abdomen, followed by 2D Cine images and prospectively ECG-gated 2D Phase contrast images prescribed at four locations along the aorta and one location in the left common carotid artery, as shown in figure (9). Imaging parameters are provided in table (7). The vessel area at each location was quantified by segmenting the 2D Cine images Yushkevich et al. (2006). The velocity at each location was quantified by manual segmentation of the phase difference images (ImageJ 1.48v;imagej.nih.gov/ij)1.48v;imagej.nih.gov/ij), extracting the mean intensity (mi) inside the contour at each time frame, and computing the following equation: u = (2048-mi)/2048*venc, where venc=velocity encoding parameter. The venc values for aorta1, 2, 3, 4, and the carotid artery were 200, 150, 150, 150, and 150 cm/s, respectively. The arc lengths between consecutive planes were computed by segmentation of the bSSFP localizer images. The centerline extraction, and path length parameterization was acquired by using Seg3D 2.3.0 sci.utah.edu/cibc-software/seg3d.html CIBC (2016) and VMTK vmtk.org Antiga et al. (2008). Finally, we utilized a periodic kernel Gaussian Process regression scheme Rasmussen (2004) to smooth out the input waveforms before we provide them as training data to the physics-informed neural networks.

2D Cine 2D Phase Contrast
Flip angle (degrees) 55 25
Repetition Time/Echo Time (ms) 51.74/1.7 10.35/6.92
Bandwidth (Hz/pixel) 765 390
Voxel size (mm3) 0.65 x 0.65 x 8 0.49 x 0.49 x 5
Temporal Resolution (ms) 29.4-34.8 20.7
Cardiac Phases 30 (retrospectively-gated) 37-43 (prospectively-gated)
Table 7: Flow through the aorta/carotid bifurcation of a healthy human subject: MRI parameters for 2D Cine images and 2D Phase contrast images used in the clinical data acquisition process.

3.3.3 Numerical results

For this case we will compare the clinically measured data and the model prediction for the Aorta point. The numerical results are presented in figure (10). Comparing the clinical measurements of area and velocity at Aorta to our model predictions, we show that there is good agreement. Further more, the predicted pressure at this point is within a range consistent with the values reported in literature for a healthy patient (i.e., ranging between 80-120mmHg). We present in figure (11) each component of the loss function, equation (9), as a functions of the number of training iterations, to exhibit the training efficiency of the method.

Figure 10: Flow through the aorta/carotid bifurcation of a healthy human subject: (a) Comparison of the neural network model predictions against the clinical cross-sectional area measurements at the test point Aorta . (b) Comparison of neural network model predictions against the clinical cross-sectional velocity measurements at the test point Aorta . (c) Predicted pressure wave at the test point Aorta .
Figure 11: Flow through the aorta/carotid bifurcation of a healthy human subject: Values of the loss functions versus number of stochastic gradient descent iterations during the training of the physics-informed neural networks. The magenda colored line denotes the continuity loss (see equations (16), (17) and (18) ). The blue colored line denotes the reconstruction loss of the area measurements. The red colored line denotes the reconstruction loss of velocity measurements. The black colored line denotes the residual loss (see equation (15)). The green colored line denotes the pressure loss, meaning how close is the pressure predicted by the network with the one produced by inserting the predicted area to the pressure-area relation (see equation (2)).

As discussed in section (2.8), we can also calibrate the and parameters of three-element Windkessel models for each of the outlets using the predicted pressure and flow rate at the carotid and Aorta points. For more details we refer the reader to table (8).

vesselParam P waveforms
Carotid 2.09e+09 2.76e-10 3.98e-02
Aorta 1.48e+08 9.23e-09 7.35e-02
Table 8: Flow through the aorta/carotid bifurcation of a healthy human subject: Windkessel model parameters discovered by adaptive grid search, and the relative error in the pressure wave predicted by the neural networks versus the wave computed by evaluating the 3-element Windkessel model using the identified parameters.

3.3.4 Discontinuous Galerkin simulation using the identified Windkessel model parameters

Using the discovered three-element Windkessel model parameters we can employ a conventional Discontinuous Galerkin solver to infer the velocity within the arterial network, and compare these with results against the reference measurements and the neural network model predictions. For the DG simulation, we choose the blood density to be equal to Kg/m and the viscosity mPas Fossan et al. (2018). Moreover, we provide the DG solver with the Windkessel and structural parameters introduced in table (9). The results are presented in figure (12).

Arterial Length Peripheral Peripheral Equilibrium cross-
segment (cm) Resistance compliance sectional area
(10 Pa s m) (10 m Pa) (Pa/m) (m)
1 4.964 - - -6.47581e+06*x + 2.47267e+06 6.91e-04*x + 2.29e-04
2 5.32 0.2285 2.759 1.37380e+08*x + 2.15121e+06 -4.46e-03*x + 2.64e-04
3 8.866 - - 7.32369e+06*x + 2.15121e+06 -5.18e-04*x + 2.64e-04
4 3.2259 0.01667 92.32 -8.42729e+06*x + 2.80053e+06 7.26e-04*x + 2.18e-04
Table 9: Flow through the aorta/carotid bifurcation of a healthy human subject: Physiological data used as simulation parameters for the Discontinuous Galerkin method for the aorta/carotid bifurcation case. Based on the 4D flow MRI measurements, we have assumed a linear tapering across the vessels’ length (i.e. the equilibrium cross-sectional area and the parameter are linear functions of the local spatial coordinate of each artery).

In section (2.1) we argued that the presented one dimensional reduced order model Sherwin et al. (2003b) constitutes an accurate approximation of the underlying fluid dynamics, but here we do actually observe some discrepancy between the model predictions and the clinically acquired data, see figure (12) for Aorta and the carotid outlets. The accuracy of the one-dimensional model has been previously validated by Reymond et. al. for a both nominal as well as patient-specific geometries Reymond et al. (2009, 2011). In both cases the authors considered arterial trees consisting of many systemic arteries including the Left Common Carotid Artery and the Thoracic Aorta. In the case of patient specific validation, they employed an average of 4 measurements of the local diameter at each cross-section and also averaged over 5-15 cardiac cycles to obtain a smoothed flow waveform. Moreover, for modeling the pulse wave propagation they utilized a form of incompressible Navier-Stokes equations that takes into consideration the wall viscous effects Reymond et al. (2009), as well as a non-linear viscoelastic pressure relation for which they tuned the viscoelastic parameters based on values appearing in the literature. Even in this more detailed case, where more complex models that can better capture the complexity of the real world are utilized, there existed some discrepancy between the predicted values and the measured ones Reymond et al. (2009, 2011). In our case, we do not pre-process the data in such a detailed manner, nor we employ such an expressive model that requires a large amount of information related to geometry, elastic properties etc. It can be seen from figure (12) that the Discontinuous Galerkin simulation that utilizes the Windkessel parameters discovered by the procedure described in section (2.8), can capture the magnitude and the wave peak timing in a favorable manner. For the case of the Left Common Carotid Artery we detect a larger amount of discrepancy and a time shift, but this maybe attributed to back-propagating elastic waves that the model in this form can not capture. Overall, considering the simplicity of the one-dimensional pulsatile flow model and the amount of provided information, it still produces favorable results. The prediction of the neural network on the other hand is very accurate at the points that measurements are provided, but slightly over-estimates the velocity at the point that the model is not trained. This discrepancy can be attributed to the simplicity of the underlying physics-informed constraints, as well as the noise corruption in the velocity and cross-sectional area measurements used to train the neural network model. Overall, for the limited amount of measurements and the noise level they contain, the performance of the deep learning model is still quite reasonable.

Figure 12: Flow through the aorta/carotid bifurcation of a healthy human subject: Comparison of the clinically acquired waveforms of blood velocity of velocity versus the predictions of the proposed physics-informed neural networks, and of a conventional Discontinuous Galerkin solver with with a three-element Windessel model outflow condition using the and parameters identified by post processing the neural network outputs. (a) Aorta . (b) Aorta . (c) Aorta . (d) Left Common Carotid.

4 Discussion

Advances in physics-informed machine learning provide the connecting link for integrating theoretical models and real-world data. This work is just one example of this general paradigm applied to the modeling and simulation of cardiovascular flows. Here, we have demonstrated how one-dimensional models of blood flow model can by seamlessly synthesized with clinical data to construct deep neural networks that can predict quantities which cannot be reliably measured in a non-invasive manner (e.g., blood pressure), by complementing the governing laws of fluid flow in compliant arteries with scattered measurements obtained by medical imaging techniques. To facilitate the efficient training of such physics-informed networks, we put forth proper non-dimensionalization and normalization techniques, as well as formulated a composite training objective that enables the consistent propagation of flow information across a networks of systemic arteries. The proposed methodology yielded favorable results across a collection of synthetic and realistic examples, and showed good agreement with state-of-the-art numerical solvers using Discontinuous Galerkin discretizations. In contrast to traditional numerical solvers, our approach does not require the precise prescription of boundary conditions and can be easily adapted to new scenarios. Moreover, we showed how a simple evaluation of the trained neural networks can provide a low-cost post-processing procedure for calibrating Windkessel model parameters for cases where traditional simulations are sought.

The direct incorporation of clinical measurements is one of the main strengths of the proposed methodology, but also one of its weaknesses. Such measurements have often very coarse resolution and may be heavily corrupted by noise. This setting introduces challenges in encouraging the neural networks to simultaneously fit the data and satisfy the underlying conservation laws, as these requirements may contradict each other in the presence of noise. Although, this did not pose a significant challenge in the reported clinical case, it may become an issue for studying smaller arterial networks for which the resolution of 4D flow MRI techniques may not suffice to produce readings with low signal-to-noise ratio. To this end, the recent work of Rudy et. al. Rudy et al. (2018) may provide a promising direction for signal-noise decomposition that is well aligned with the physics-informed deep learning framework presented in this work.

In terms of future work, we are mainly focused on testing and validating the proposed method on more clinical cases. To this end, we consider extensions to larger arterial networks having more complex geometry, e.g., many bifurcation and junction points and an increased number of vessels. Furthermore, we plan to test the effectiveness of neural network architectures, including self-supervised learning with convolutional networks

Zhu et al. (2019), in order to decrease the number of parameters, and, consequently, the wall clock time of the analysis. Another direction for future development is related to extending our methods for handling more realistic arterial wall mechanics and blood rheology models, accounting for visco-elastic and non-Newtonian effects.

Acknowledgments

The authors would like to thank Veronica Aramendia-Vidaurreta and Brianna Moon for assisting in the acquisition of carotid/aorta MRI data. G.K., Y.Y. and P.P. would like to acknowledge support from the US Department of Energy under the Advanced Scientific Computing Research program (grant DE-SC0019116), the Defense Advanced Research Projects Agency under the Physics of Artificial Intelligence program (grant HR00111890034). E.H. received support from the Eunice Kennedy Shriver National Institute of Child Health and Human Development (U01-HD087180), National Institute of Biomedical Imaging and Bio-engineering (T32-EB009384) and National Science Foundation (DGE-1321851).

References

  • Aaslid et al. (1982) Aaslid, R., Markwalder, T.-M., Nornes, H., 1982. Noninvasive transcranial doppler ultrasound recording of flow velocity in basal cerebral arteries. Journal of neurosurgery 57 (6), 769–774.
  • Abadi et al. (2016) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al., 2016. Tensorflow: a system for large-scale machine learning. In: OSDI. Vol. 16. pp. 265–283.
  • Alastruey et al. (2008) Alastruey, J., Moore, S., Parker, K., David, T., Peiró, J., Sherwin, S., 2008. Reduced modelling of blood flow in the cerebral circulation: Coupling 1-d, 0-d and cerebral auto-regulation models. International journal for numerical methods in fluids 56 (8), 1061–1067.
  • Antiga et al. (2008) Antiga, L., Piccinelli, M., Botti, L., Ene-Iordache, B., Remuzzi, A., Steinman, D. A., 2008. An image-based modeling framework for patient-specific computational hemodynamics. Medical & biological engineering & computing 46 (11), 1097.
  • Barker et al. (1985) Barker, A. T., Jalinous, R., Freeston, I. L., 1985. Non-invasive magnetic stimulation of human motor cortex. The Lancet 325 (8437), 1106–1107.
  • Chan et al. (2007) Chan, C., Poon, C., Wong, R. C., Zhang, Y., 2007. A hybrid body sensor network for continuous and long-term measurement of arterial blood pressure. In: 2007 4th IEEE/EMBS International Summer School and Symposium on Medical Devices and Biosensors. IEEE, pp. 121–123.
  • Chen et al. (2013) Chen, P., Quarteroni, A., Rozza, G., 2013. Simulation-based uncertainty quantification of human arterial network hemodynamics. International journal for numerical methods in biomedical engineering 29 (6), 698–721.
  • CIBC (2016) CIBC, 2016. Seg3D: Volumetric Image Segmentation and Visualization. Scientific Computing and Imaging Institute (SCI), Download from: http://www.seg3d.org.
  • Cockburn et al. (2012) Cockburn, B., Karniadakis, G. E., Shu, C.-W., 2012. Discontinuous Galerkin methods: theory, computation and applications. Vol. 11. Springer Science & Business Media.
  • Duchi et al. (2011)

    Duchi, J., Hazan, E., Singer, Y., 2011. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research 12 (Jul), 2121–2159.

  • Edelstein et al. (1980) Edelstein, W. A., Hutchison, J. M., Johnson, G., Redpath, T., 1980. Spin warp nmr imaging and applications to human whole-body imaging. Physics in medicine & biology 25 (4), 751.
  • Famiglietti and Wood (1994) Famiglietti, J., Wood, E., 1994. Multiscale modeling of spatially variable water and energy balance processes. Water Resources Research 30 (11), 3061–3078.
  • Figueroa et al. (2006) Figueroa, C. A., Vignon-Clementel, I. E., Jansen, K. E., Hughes, T. J., Taylor, C. A., 2006. A coupled momentum method for modeling blood flow in three-dimensional deformable arteries. Computer methods in applied mechanics and engineering 195 (41-43), 5685–5706.
  • Formaggia et al. (2010) Formaggia, L., Quarteroni, A., Veneziani, A., 2010. Cardiovascular Mathematics: Modeling and simulation of the circulatory system. Vol. 1. Springer Science & Business Media.
  • Fossan et al. (2018) Fossan, F. E., Mariscal-Harana, J., Alastruey, J., Hellevik, L. R., 2018. Optimization of topological complexity for one-dimensional arterial blood flow models. Journal of the Royal Society Interface 15 (149), 20180546.
  • Glorot and Bengio (2010) Glorot, X., Bengio, Y., 2010. Understanding the difficulty of training deep feedforward neural networks. In: Proceedings of the thirteenth international conference on artificial intelligence and statistics. pp. 249–256.
  • Grinberg et al. (2011) Grinberg, L., Cheever, E., Anor, T., Madsen, J. R., Karniadakis, G., 2011. Modeling blood flow circulation in intracranial arterial networks: a comparative 3d/1d simulation study. Annals of biomedical engineering 39 (1), 297–309.
  • Grinberg and Karniadakis (2008) Grinberg, L., Karniadakis, G. E., 2008. Outflow boundary conditions for arterial networks with multiple outlets. Annals of biomedical engineering 36 (9), 1496–1514.
  • Haris et al. (2014) Haris, M., Singh, A., Cai, K., Kogan, F., McGarvey, J., DeBrosse, C., Zsido, G. A., Witschey, W. R., Koomalsingh, K., Pilla, J. J., et al., 2014. A technique for in vivo mapping of myocardial creatine kinase metabolism. Nature medicine 20 (2), 209.
  • Ismail et al. (2013) Ismail, M., Wall, W. A., Gee, M. W., 2013. Adjoint-based inverse analysis of windkessel parameters for patient-specific vascular models. Journal of Computational Physics 244, 113–130.
  • Kak et al. (2002) Kak, A. C., Slaney, M., Wang, G., 2002. Principles of computerized tomographic imaging. Medical Physics 29 (1), 107–107.
  • Kett-White et al. (2002) Kett-White, R., Hutchinson, P. J., Al-Rawi, P. G., Gupta, A. K., Pickard, J. D., Kirkpatrick, P. J., 2002. Adverse cerebral events detected after subarachnoid hemorrhage using brain oxygen and microdialysis probes. Neurosurgery 50 (6), 1213–1222.
  • Kingma and Ba (2014) Kingma, D. P., Ba, J., 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  • Ma et al. (2018) Ma, Y., Choi, J., Hourlier-Fargette, A., Xue, Y., Chung, H. U., Lee, J. Y., Wang, X., Xie, Z., Kang, D., Wang, H., et al., 2018. Relation between blood pressure and pulse wave velocity for human arteries. Proceedings of the National Academy of Sciences 115 (44), 11144–11149.
  • Markl et al. (2012) Markl, M., Frydrychowicz, A., Kozerke, S., Hope, M., Wieben, O., 2012. 4d flow mri. Journal of Magnetic Resonance Imaging 36 (5), 1015–1036.
  • Matthys et al. (2007) Matthys, K. S., Alastruey, J., Peiró, J., Khir, A. W., Segers, P., Verdonck, P. R., Parker, K. H., Sherwin, S. J., 2007. Pulse wave propagation in a model human arterial network: assessment of 1-d numerical simulations against in vitro measurements. Journal of biomechanics 40 (15), 3476–3486.
  • Mitchell et al. (2011) Mitchell, G. F., van Buchem, M. A., Sigurdsson, S., Gotal, J. D., Jonsdottir, M. K., Kjartansson, Ó., Garcia, M., Aspelund, T., Harris, T. B., Gudnason, V., et al., 2011. Arterial stiffness, pressure and flow pulsatility and brain structure and function: the age, gene/environment susceptibility–reykjavik study. Brain 134 (11), 3398–3407.
  • Olufsen (1999) Olufsen, M. S., 1999. Structured tree outflow condition for blood flow in larger systemic arteries. American journal of physiology-Heart and circulatory physiology 276 (1), H257–H268.
  • O’Rourke (1995) O’Rourke, M., 1995. Mechanical principles in arterial disease. Hypertension 26 (1), 2–9.
  • Pant et al. (2014) Pant, S., Fabrèges, B., Gerbeau, J.-F., Vignon-Clementel, I., 2014. A methodological paradigm for patient-specific multi-scale cfd simulations: from clinical measurements to parameter estimates for individual analysis. International journal for numerical methods in biomedical engineering 30 (12), 1614–1648.
  • Paszke et al. (2017)

    Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A., 2017. Automatic differentiation in pytorch.

  • Perdikaris et al. (2015) Perdikaris, P., Grinberg, L., Karniadakis, G. E., 2015. An effective fractal-tree closure model for simulating blood flow in large arterial networks. Annals of biomedical engineering 43 (6), 1432–1442.
  • Perdikaris and Karniadakis (2014) Perdikaris, P., Karniadakis, G. E., 2014. Fractional-order viscoelasticity in one-dimensional blood flow models. Annals of biomedical engineering 42 (5), 1012–1023.
  • Plein et al. (2001) Plein, S., Bloomer, T. N., Ridgway, J. P., Jones, T. R., Bainbridge, G. J., Sivananthan, M. U., 2001. Steady-state free precession magnetic resonance imaging of the heart: comparison with segmented k-space gradient-echo imaging. Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine 14 (3), 230–236.
  • Raissi et al. (2019) Raissi, M., Perdikaris, P., Karniadakis, G., 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686–707.
  • Raissi et al. (2017a) Raissi, M., Perdikaris, P., Karniadakis, G. E., 2017a. Machine learning of linear differential equations using gaussian processes. Journal of Computational Physics 348, 683–693.
  • Raissi et al. (2017b) Raissi, M., Perdikaris, P., Karniadakis, G. E., 2017b. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint arXiv:1711.10561.
  • Raissi et al. (2017c) Raissi, M., Perdikaris, P., Karniadakis, G. E., 2017c. Physics informed deep learning (part ii): data-driven discovery of nonlinear partial differential equations. arXiv preprint arXiv:1711.10566.
  • Rasmussen (2004) Rasmussen, C. E., 2004. Gaussian processes in machine learning. In: Advanced lectures on machine learning. Springer, pp. 63–71.
  • Revie et al. (2013) Revie, J. A., Stevenson, D. J., Chase, J. G., Hann, C. E., Lambermont, B. C., Ghuysen, A., Kolh, P., Shaw, G. M., Heldmann, S., Desaive, T., 2013. Validation of subject-specific cardiovascular system models from porcine measurements. Computer methods and programs in biomedicine 109 (2), 197–210.
  • Reymond et al. (2011) Reymond, P., Bohraus, Y., Perren, F., Lazeyras, F., Stergiopulos, N., 2011. Validation of a patient-specific one-dimensional model of the systemic arterial tree. American Journal of Physiology-Heart and Circulatory Physiology 301 (3), H1173–H1182.
  • Reymond et al. (2009) Reymond, P., Merenda, F., Perren, F., Rufenacht, D., Stergiopulos, N., 2009. Validation of a one-dimensional model of the systemic arterial tree. American Journal of Physiology-Heart and Circulatory Physiology 297 (1), H208–H222.
  • Rose et al. (1982) Rose, G. A., Blackburn, H., Gillum, R., Prineas, R., et al., 1982. Cardiovascular survey methods. Vol. 56. Geneva, Switzerland; WHO.
  • Ruder (2016) Ruder, S., 2016. An overview of gradient descent optimization algorithms. arXiv preprint arXiv:1609.04747.
  • Rudy et al. (2018) Rudy, S. H., Kutz, J. N., Brunton, S. L., 2018. Deep learning of dynamics and signal-noise decomposition with time-stepping constraints. arXiv preprint arXiv:1808.02578.
  • Schiavazzi et al. (2017) Schiavazzi, D. E., Baretta, A., Pennati, G., Hsia, T.-Y., Marsden, A. L., 2017. Patient-specific parameter estimation in single-ventricle lumped circulation models under uncertainty. International journal for numerical methods in biomedical engineering 33 (3), e02799.
  • Sherwin et al. (2003a) Sherwin, S., Formaggia, L., Peiro, J., Franke, V., 2003a. 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), 673–700.
  • Sherwin et al. (2003b) Sherwin, S., Franke, V., Peiró, J., Parker, K., 2003b. One-dimensional modelling of a vascular network in space-time variables. Journal of Engineering Mathematics 47 (3-4), 217–250.
  • Spilker et al. (2007) Spilker, R. L., Feinstein, J. A., Parker, D. W., Reddy, V. M., Taylor, C. A., 2007. Morphometry-based impedance boundary conditions for patient-specific modeling of blood flow in pulmonary arteries. Annals of biomedical engineering 35 (4), 546–559.
  • Spilker and Taylor (2010) Spilker, R. L., Taylor, C. A., 2010. Tuning multidomain hemodynamic simulations to match physiological measurements. Annals of biomedical engineering 38 (8), 2635–2648.
  • Stefanovska (1999) Stefanovska, A., 1999. Physics of the human cardiovascular system. Contemporary Physics 40 (1), 31–55.
  • Stein (1987) Stein, M., 1987. Large sample properties of simulations using latin hypercube sampling. Technometrics 29 (2), 143–151.
  • Tartakovsky et al. (2018) Tartakovsky, A. M., Marrero, C. O., Tartakovsky, D., Barajas-Solano, D., 2018. Learning parameters and constitutive relationships with physics informed deep neural networks. arXiv preprint arXiv:1808.03398.
  • Urquiza et al. (2006) Urquiza, S., Blanco, P., Vénere, M., Feijóo, R., 2006. Multidimensional modelling for the carotid artery blood flow. Computer Methods in Applied Mechanics and Engineering 195 (33-36), 4002–4017.
  • Westerhof et al. (2009) Westerhof, N., Lankhaar, J.-W., Westerhof, B. E., Feb 2009. The arterial windkessel. Medical & Biological Engineering & Computing 47 (2), 131–141.
    URL https://doi.org/10.1007/s11517-008-0359-2
  • Yang and Perdikaris (2018) Yang, Y., Perdikaris, P., 2018. Adversarial uncertainty quantification in physics-informed neural networks. arXiv preprint arXiv:1811.04026.
  • Yu et al. (1998) Yu, Y.-C., Boston, J. R., Simaan, M. A., Antaki, J. F., 1998. Estimation of systemic vascular bed parameters for artificial heart control. IEEE Transactions on Automatic Control 43 (6), 765–778.
  • Yushkevich et al. (2006) Yushkevich, P. A., Piven, J., Hazlett, H. C., Smith, R. G., Ho, S., Gee, J. C., Gerig, G., 2006. User-guided 3d active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage 31 (3), 1116–1128.
  • Zhu and Zabaras (2018) Zhu, Y., Zabaras, N., 2018. Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics 366, 415–447.
  • Zhu et al. (2019) Zhu, Y., Zabaras, N., Koutsourelakis, P.-S., Perdikaris, P., 2019. Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data. arXiv preprint arXiv:1901.06314.