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), biomedical 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), biophysics 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 invivo 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 KettWhite et al. (2002). These limitations motivate the use of noninvasive measurement techniques such as biomedical 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 threedimensional 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 crosssectional 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 noninvasive 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 invivo and noninvasive manner have motivated the use of computer simulations and computational fluid dynamics models to predict them insilico. Advances in algorithms and computing resources now allow us to perform detailed flow simulations in complex patientspecific arterial topologies using threedimensional (3D) and/or onedimensional (1D) formulations of the unsteady NavierStokes 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 invitro experiments Matthys et al. (2007) as well as invivo 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 noninvasive 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 physicsinformed 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 noninvasive 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 crosssections), (ii) satisfy the underlying physical conservation laws as described by a reduced onedimensional 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, pretrained 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 postprocessing 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 onedimensional 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 invivo 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
Pulse wave propagation in arterial networks can be effectively modeled using onedimensional (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 crosssection. Following Sherwin et al. (2003a); Formaggia et al. (2010) we consider a reduced form of the incompressible NavierStokes 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 crosssectional 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 invitro and invivo 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 crosssectional 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 Physicsinformed 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 noninvasive 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 noninvasively 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.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 latinhypercube 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
Onedimensional 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.
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 crosssectional area. The second summation encourages the model to satisfy the continuity equation on the bifurcation point. The third summation represents the physicsinformed 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 physicsinformed 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 Nondimensionalization and normalization
In equation (1), the order of magnitude of the different physical quantities, velocity, crosssectional 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 nondimensionalization 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 nondimensionalization we introduce the characteristic variables, which are commonly used in multiscale 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 crosssectional 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 nondimensional and normalized form, all the variables and inputs are scaled to order . This is what effectively enables the training of our physicsinformed 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 nondimensional continuity conditions, by inserting the nondimensionalizing 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 nondimensional 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 threeelement 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 threeelement Windkessel parameters for a coupled 3D/0D system by implementing an inverse analysis based on adjoint methods, using the patientspecific pressure information as a target. Spilker et. al. Spilker and Taylor (2010) proposed a multidomain method that gets coupled with a reduced order, a threeelement 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. YihChoung 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 zerodimensional model with known pressure, at specific points, to recover the values of threeelemnent 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 physicsinformed neural networks (as described in section (2)) does not require the apriori specification of Windkesseltype 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 postprocessing 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 physicsinformed 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 Yshaped 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 timeseries 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 inhouse 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.36e05 
2  0.7  0.5251  0.3428  5.42e+08  1.81e06 
3  0.67  0.2702  0.6661  6.96e+07  1.36e05 
3.1.3 Neural network approximation setup
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 tradeoff of computational cost and predictive accuracy. The neural networks take the scaled inputs and predict the nondimensional outputs , as discussed in section (2.7). We are providing each neural network with velocity and crosssectional 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 finetuning the neural network parameters.3.1.4 Comparison of model predictions and reference solution
We first demonstrate the importance of the nondimensionalization 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 nondimensionalization tools significantly improve the model training capabilities and, thus, provide accurate predictions.
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 physicinformed neural networks model. The prediction error is further quantified in the systematic studies presented in section (3.1.5).
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.
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 hyperparameter settings. All numerical experiments are performed for the canonical Yshaped 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.42e02  4.99e02  5.65e02  3.92e02  1.07e01  5.52e02  5.34e02  3.54e02  6.41e02  4.57e02 
4.62e02  4.04e02  2.94e02  3.78e02  3.48e02  1.32e01  6.01e02  3.00e02  4.78e02  4.63e02 
3.68e02  4.48e02  4.37e02  5.81e02  1.23e01  5.80e02  4.06e02  2.33e02  3.02e02  3.96e02 
5.70e02  3.61e02  3.85e02  3.83e02  4.82e02  4.85e02  5.46e02  6.01e02  3.84e02  7.43e02 
5.82e02  5.09e02  4.13e02  6.66e02  1.31e01  7.31e02  6.09e02  5.59e02  3.50e02  3.47e02 
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 noisefree 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 physicsinformed constraint on the PDEs residual can effectively regularize the training process and safeguard against overfitting. 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.14e01  8.41e01  7.04e01  1.07e+00 
3  4.09e01  3.08e01  3.68e01  2.39e01 
5  5.18e01  2.35e01  1.03e01  2.49e01 
7  4.42e01  3.88e02  3.73e02  3.56e02 
3.1.6 Windkessel parameter identification by adaptive grid search
As described in section (2.8), our method provides an effective mechanism for c//;/alibrating Windkessel model parameters as a simple postprocessing 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.89e11  3.36e02  
1.84e+09  7.36e11  3.67e02 
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 crosssectional 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 crosssectional area and velocity using an inhouse 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.14e05 
2  6.66638      2.60e+07  2.21e05 
3  6.99352      2.63e+07  2.17e05 
4  14.7735  0.3133  16.62  2.82e+07  1.97e05 
5  14.9503  0.1654  31.49  2.71e+07  2.08e05 
6  13.6421  0.1682  30.96  2.67e+07  2.12e05 
7  13.4384  0.2086  2.092  2.87e+07  1.92e05 
3.2.3 Neural network approximation setup
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 nondimensional 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 Latinhypercube 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 threeelement 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.
vesselParam  P waveforms  

1.36e+09  3.78e09  7.92e03  
1.64e+09  2.22e10  3.29e02  
1.22e+09  2.96e09  8.26e03  
2.69e+09  1.84e09  2.56e02 
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 crosssectional 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).
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 headfirst, supine on the table. For the 3lead 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 SteadyState Free Precession (bSSFP) localizer of the neck through abdomen, followed by 2D Cine images and prospectively ECGgated 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 = (2048mi)/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/cibcsoftware/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 physicsinformed 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.434.8  20.7 
Cardiac Phases  30 (retrospectivelygated)  3743 (prospectivelygated) 
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 80120mmHg). 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.
As discussed in section (2.8), we can also calibrate the and parameters of threeelement 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.76e10  3.98e02 
Aorta  1.48e+08  9.23e09  7.35e02 
3.3.4 Discontinuous Galerkin simulation using the identified Windkessel model parameters
Using the discovered threeelement 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.91e04*x + 2.29e04 
2  5.32  0.2285  2.759  1.37380e+08*x + 2.15121e+06  4.46e03*x + 2.64e04 
3  8.866      7.32369e+06*x + 2.15121e+06  5.18e04*x + 2.64e04 
4  3.2259  0.01667  92.32  8.42729e+06*x + 2.80053e+06  7.26e04*x + 2.18e04 
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 onedimensional model has been previously validated by Reymond et. al. for a both nominal as well as patientspecific 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 crosssection and also averaged over 515 cardiac cycles to obtain a smoothed flow waveform. Moreover, for modeling the pulse wave propagation they utilized a form of incompressible NavierStokes equations that takes into consideration the wall viscous effects Reymond et al. (2009), as well as a nonlinear 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 preprocess 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 backpropagating elastic waves that the model in this form can not capture. Overall, considering the simplicity of the onedimensional 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 overestimates the velocity at the point that the model is not trained. This discrepancy can be attributed to the simplicity of the underlying physicsinformed constraints, as well as the noise corruption in the velocity and crosssectional 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.
4 Discussion
Advances in physicsinformed machine learning provide the connecting link for integrating theoretical models and realworld 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 onedimensional 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 noninvasive 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 physicsinformed networks, we put forth proper nondimensionalization 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 stateoftheart 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 lowcost postprocessing 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 signaltonoise ratio. To this end, the recent work of Rudy et. al. Rudy et al. (2018) may provide a promising direction for signalnoise decomposition that is well aligned with the physicsinformed 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 selfsupervised 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 viscoelastic and nonNewtonian effects.Acknowledgments
The authors would like to thank Veronica AramendiaVidaurreta 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 DESC0019116), 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 (U01HD087180), National Institute of Biomedical Imaging and Bioengineering (T32EB009384) and National Science Foundation (DGE1321851).
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 largescale 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 1d, 0d and cerebral autoregulation models. International journal for numerical methods in fluids 56 (8), 1061–1067.
 Antiga et al. (2008) Antiga, L., Piccinelli, M., Botti, L., EneIordache, B., Remuzzi, A., Steinman, D. A., 2008. An imagebased modeling framework for patientspecific computational hemodynamics. Medical & biological engineering & computing 46 (11), 1097.
 Barker et al. (1985) Barker, A. T., Jalinous, R., Freeston, I. L., 1985. Noninvasive 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 longterm 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. Simulationbased 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 wholebody 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., VignonClementel, I. E., Jansen, K. E., Hughes, T. J., Taylor, C. A., 2006. A coupled momentum method for modeling blood flow in threedimensional deformable arteries. Computer methods in applied mechanics and engineering 195 (4143), 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., MariscalHarana, J., Alastruey, J., Hellevik, L. R., 2018. Optimization of topological complexity for onedimensional 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. Adjointbased inverse analysis of windkessel parameters for patientspecific 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.
 KettWhite et al. (2002) KettWhite, R., Hutchinson, P. J., AlRawi, 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., HourlierFargette, 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 1d 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 physiologyHeart 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., VignonClementel, I., 2014. A methodological paradigm for patientspecific multiscale 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 fractaltree 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. Fractionalorder viscoelasticity in onedimensional 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. Steadystate free precession magnetic resonance imaging of the heart: comparison with segmented kspace gradientecho 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. Physicsinformed 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): Datadriven 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): datadriven 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 subjectspecific 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 patientspecific onedimensional model of the systemic arterial tree. American Journal of PhysiologyHeart 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 onedimensional model of the systemic arterial tree. American Journal of PhysiologyHeart 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 signalnoise decomposition with timestepping constraints. arXiv preprint arXiv:1808.02578.
 Schiavazzi et al. (2017) Schiavazzi, D. E., Baretta, A., Pennati, G., Hsia, T.Y., Marsden, A. L., 2017. Patientspecific parameter estimation in singleventricle 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 (67), 673–700.
 Sherwin et al. (2003b) Sherwin, S., Franke, V., Peiró, J., Parker, K., 2003b. Onedimensional modelling of a vascular network in spacetime variables. Journal of Engineering Mathematics 47 (34), 217–250.
 Spilker et al. (2007) Spilker, R. L., Feinstein, J. A., Parker, D. W., Reddy, V. M., Taylor, C. A., 2007. Morphometrybased impedance boundary conditions for patientspecific 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., BarajasSolano, 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 (3336), 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/s1151700803592  Yang and Perdikaris (2018) Yang, Y., Perdikaris, P., 2018. Adversarial uncertainty quantification in physicsinformed 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. Userguided 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. Physicsconstrained deep learning for highdimensional surrogate modeling and uncertainty quantification without labeled data. arXiv preprint arXiv:1901.06314.
Comments
There are no comments yet.