Verification, validation and uncertainty quantification in cardiovascular CFD

Introduction
Image-based CFD of cerebral aneurysms is subject to appreciable uncertainties, notably the flow rate boundary conditions, which are rarely available clinically. This is exacerbated by recent reports of prevalent high-frequency flow instabilities [1], which are more sensitive to inflow conditions [2]. Comprehensive uncertainty quantification (UQ) is not practical using such high-fidelity CFD, so here we aimed to demonstrate that coarsened CFD could reasonably capture the flow instabilities, and thus allow “accurate UQ using inaccurate computational models” [3]

Methods
14 bifurcation aneurysm cases from the Aneurisk database were meshed to generate, on average, 3 million (3M) tetrahedral elements. Pulsatile simulations were performed using second-order velocity shape functions (so equivalent to ~24M linear tetrahedra) with 20,000 time-steps/cycle, using a representative waveform with cycle-averaged inflow rate scaled to inlet diameter [4]. For one case with strong flow instabilities, simulations were also run with cycle-averaged inflow varied by a factor of 8.

Coarser meshes were generated by increasing the nominal tetrahedral side length by factors of 1.5 and 2.0, resulting in meshes averaging 1M and 0.5M (linear) tetrahedra, respectively. These were run at the same flow rates, with 5,000 steps/cycle, requiring 10 and 4 hrs/cycle vs. 150 for the fine simulations. For all cases we computed sac-averaged values of: time-averaged wall shear stress, normalized to the parent artery (WSS*); oscillatory shear index (OSI); and spectral power index (SPI), a recently-introduced metric for quantifying high frequency flow instabilities [2].

Results
Per the top row of the figure below, coarse results were significantly correlated with fine ones. For SPI, slopes were significantly<1, confirming a consistent (and expected) underestimation of SPI by coarse CFD. These trends are reflected in the bottom row of the figure, where fine and coarse CFD show the same variations with flow rate.



Discussion
Multifidelity UQ allows coarse simulations to deviate from fine ones, as long as “they exhibit some sort of dependence” [3].  Consistent with an earlier study focused on only 3 aneurysm cases [5], the present study confirms this dependency, even for a sensitive quantity like SPI.  This encourages us that comprehensive UQ can be performed effectively, since each fine simulation can be replaced by at least 40 coarse ones. Further coarsening of elements caused some failed meshes; however, reductions might still be achievable by further coarsening and/or truncation of inflow/outflow branches.

Acknowledgements 
Heart & Stroke Foundation grant G-16-00012564

References
1.  Valen-Sendstad K, et al. AJNR Am J Neuroradiol 2014;35(3):536-43.
2.  Khan MO, et al. J Biomech 2017;52:179-182.
3.  Koutsourelakis PS. SIAM J Sci Comput 2009;31(5):3274-3300.
4.  Valen-Sendstad K, et al. Ann Biomed Eng 2015;43(6):1422-31.
5.  Khan MO, et al. AJNR Am J Neuroradiol 2015;36(7):1310-6.

Introduction

Doppler ultrasound (DUS) is a powerful tool for detecting life risking pathologies during pregnancy by providing real time flow field images of blood vessels in the placenta. However, a professional reading of the DUS in terms of blood flow field is a necessity, which otherwise pathologies might be overlooked or mistaken. Therefore, the objective of the presented study is to explore the flow field of different umbilical cord (UC) models, which represent healthy and pathological conditions, using DUS, compare the results with Computational Fluid Dynamics (CFD) data simulated on the same models, and to examine the correlation between the two methods in order to provide an assessment of the ability of the DUS to identify defects and complications associated with UC flow.

Methods

Three different approaches were applied in order to compare their flow field outcome using the same UC models, two measurements methods and one computational. The measurements obtained within a dynamic simulator [1], which producing a controlled pulsatile flow (Fig. 1A) into a pre-determined vessel model: using (1)  a built-in flow sensor, (2) a portable DUS machine, while the tube is located in a water container that mimics a "uterus". The CFD (3) Ansys Fluent, was used to simulate the same carried out experiments for outlining the transient flow field pattern inside the UC using three-dimensional (3D) models. Three separate vessel models were tested: straight, helical and knotted tubes.

Results

Two types of correlations were performed: Velocity profile and contours obtained from CFD and DUS Measurements. The straight tube correlated in the range of 5-25%, the helical tube, showed correlation with 2-16% variation. CFD contours observed at particular time steps, the captured DUS image represent an unknown time step (Fig 1). The knotted tube showed correlation range from 7%-60% variation.

Discussion

The obtained results showed correlation for the straight and helical models, the knotted tube, which has a complex geometric structure showed poor correlation. Velocity transient profile of the straight and helical tubes were able to correlate with an accepted error difference, the knotted tube gave inconsistent correlation for all points of interest, with no apparent pattern. For all models, velocity contour images yield poor correlations, as the DUS provided limited color map range.

It was concluded that pathologies with complex geometric structure may misinterpreted, as it requires skillful handling of the device. Pathologies with complex geometric structure are even more quality sensitive and may misinterpreted. Thus, it is recommended to proceed further with this research in order to establish correlation standards and to draw outline of potential limitation in DUS flow field interpretation.

References

1. Agranati O. & Naftali S. "Dynamic simulator for umbilical flow assessment using Doppler Ultrasound", Final project book, 2015, Afeka College.

Introduction

Aortic valve disease (AVD) is a condition in which the aortic valve displays reduced function due to damage or disease. While pulsatile blood flow in a healthy aorta can typically be considered laminar, evidence suggests that in the presence of valve pathology, aortic flow can become highly disturbed and turbulent near peak systolic flow [1,2].

There have been numerous computational fluid dynamic (CFD) studies of aortic haemodynamics, however, the majority of previous studies have made the assumption of laminar flow. Of the numerical models which consider turbulence, most assume fully turbulent flow with only a few exceptions [3,4]. To better understand the role of turbulence in the initiation and progression of aortic disease, laminar to turbulence transition must be considered. Therefore, the aim of this study is to simulate transitional flow in a human aorta with severe aortic stenosis using a scale-resolved turbulence model, large-eddy simulation (LES). Furthermore, results obtained from Reynolds-averaged Navier-Stokes (RANS)-based models will be used for comparison.

Methods

A patient with severe aortic stenosis was identified as suitable for the project, showing Reynolds numbers in the ascending aorta that are associated with turbulence. Patient-specific aortic geometry was reconstructed from magnetic resonance imaging (MRI) data. The model includes the ascending aorta, aortic arch, supra-aortic branches and the proximal descending and thoracic aorta.

Simulations were carried out using the open source software, OpenFOAM. LES was used as this method is capable of capturing transitional and turbulent flow features. RANS-based simulations employing the shear stress transport (SST) transitional model were also performed, together with a laminar flow simulation. Blood was assumed to be incompressible and Newtonian, with density 1060 kg/m³ and dynamic viscosity 3.5·10^-3 Pa·s. Phase contrast MR images were processed with an in-house Matlab code to generate patient-specific 3D pulsatile velocity profiles imposed at the model inlet. A three-element Windkessel model was applied at each of the four model outlets.

Results

The methodology has been implemented within OpenFOAM, providing initial results. Differences were observed between the results of LES, RANS-based SST transitional model and laminar flow simulations, highlighting the importance in selecting an appropriate flow model for aortic flow in the presence of valve pathology.

Discussion

Detailed analysis and comparisons are currently ongoing and will be presented at the conference. Due to the highly transient nature of turbulence in the aorta, it is critical to consider the temporal development of turbulence and relaminarization during a cardiac cycle.

References

1.           Stein, P.D., and Sabbah, H.N., (1976). Circ Res, 39(1) p58-65

2.           Tan, F.P.P., et al., (2011). Cardiovasc Eng Technol, 3(1) p123-135

3.           Lantz, J., et al., (2013). J Biomech, 46(11) p1851-1858

4.           Kousera, C.A., et al., (2012). J Biomech Eng, 135(1) p11003

Introduction

            In recent years, CFD has evolved from a research tool into a promising tool for medical device evaluation and clinical diagnosis. In this context, it is critical to validate CFD solvers prior to their application. However, guidelines for CFD solver validation in biomedical applications do not exist. To this end, we present a validation methodology that systematically quantifies the accuracy of CFD solvers in their intended biomedical use.

Methods

            Our validation methodology is based on the concepts of the ASME Verification & Validation 20 Standard.[1] The goal is to quantify the CFD model error (δ_model), an estimate of “how good” the assumptions of the solver are at simulating the physics of the flow. As shown in Fig. 1A, sources of error in both simulation and experiment are identified and parsed out to quantify δ_model as:

δ_model = E ± u_val                 

where error (E) is:

E = S  ̶  D                           

the difference between simulation (S) and experimental (D) results, and model uncertainty (u_val) is:

sqrt(u_num^2 + u_input^2 + u_D^2)

RMS of numerical uncertainty (u_num), input parameter uncertainty (u_input) and experimental uncertainty (u_D). u_num is the uncertainty due to spatial and temporal discretizations, u_input is the uncertainty due to errors in simulation input parameters, and u_D is the uncertainty in experimental measurements.

            The validation methodology was applied to quantify the accuracy of STAR-CCM+ (CD-adapco, Melville, NY) in a clinically-relevant context by simulating blood flow (density: 1056 kg/m³, viscosity: 0.00415 Pa-s) in a 9-mm patient-specific intracranial aneurysm located at the internal carotid artery under pulsatile flow conditions (flow rate: 256 ml/min). 2D particle image velocimetry (PIV) in an identical silicone model provided experimental velocity measurements. Velocity magnitude from simulation (CFD) and PIV results were used to calculated E, u_val and δ_model.

Results

            Qualitative comparison (Fig. 1B) showed good match between velocity fields in CFD and PIV. The accuracy of STAR-CCM+ was shown by plotting δ_model along the intersecting line (Fig. 1C). The highest magnitude of δ_model occurred at x=0.89 with large negative E (-4.63 cm/s) and small u_val (0.79 cm/s), meaning that at this location, CFD under-predicted velocity as compared to the experiment. Average δ_model in STAR-CCM+ was 5.63%±5.49% along the intersecting line (Fig. 1D).

Figure1: Validation analysis and results.

Discussion

            Our validation methodology offers a streamlined workflow to quantify the accuracy of CFD solvers. δ_model provides the range of error of a CFD solver in simulating flow in its context of use. This range could be crucial while choosing the best solver among candidates for their proposed clinical application.

Acknowledgements

            NIH (R01-NS-091075), Toshiba Medical System Corporation.

References

1. ASME V&V 20 (2009)

Introduction: Patient-specific computational fluid dynamics (CFD) simulations are progressively becoming a powerful tool in surgical planning for the Fontan procedure, a palliative heart surgery for single-ventricle patients resulting in a total cavopulmonary connection (TCPC). In these simulations, parabolic profiles are commonly used for velocity boundary conditions (BCs). One major drawback of parabolic velocity profiles is the lack of consideration of flow pulsatility. Therefore, more sophisticated numerical models, such as the Womersley model, have been developed to account for this shortcoming (1). However, previous literature has shown that the Womersley model is not superior to the parabolic model in comparison with the hemodynamic assessment by using the real velocity profile, despite the fact that the former requires more programming and computational cost than the latter (2). This study investigated the effectiveness of the Womersley model in simulating Fontan hemodynamics.

Method: Simulations were conducted for ten Fontan patients with four different inlet velocity profiles: (1) parabolic velocity profile,( 2) Womersley velocity profile, (3) parabolic velocity profile with inlet extensions (to develop flow before entering the TCPC), and (4) the real, patient-specific velocity profile. Each simulation was regulated such that each vessel’s flow waveform was identical to the patient-specific flow rate obtained from phase-contrast MRI and that they varied only in terms of their spatial velocity profiles at the inlets. We quantified the hemodynamic results and discrepancies of the simulations by examining three main metrics: wall shear stress, power loss, and hepatic fluid distribution, all of which have shown potential relevance in Fontan surgical planning.

Results: The parabolic and Womersley simulations resulted in a similar level of difference in all three metrics when compared to the simulation employing the patient-specific velocity profile. An example is shown in Fig. 1.

Both the parabolic and Womersley simulation results showed large discrepancies (4-20% depending on the metric) when compared to the real, patient-specific velocity profile simulation results. The simulations employing extensions at the inlet vessels with parabolic BCs were seen to most closely approximate the patient-specific profile results.

Discussion: In comparison with the real velocity profiles, simulations using parabolic and Womersley velocity profiles show significant differences in calculated Fontan hemodynamics. These results are relevant to the future of Fontan surgical planning, particularly in relation to the choice of BCs used in these simulations.

Acknowledgement: This study was supported by the National Heart, Lung, and Blood Institute Grants HL67622 and HL098252. The authors acknowledge the use of ANSYS software which was provided through an Academic Partnership between ANSYS, Inc. and the Cardiovascular Fluid Mechanics Lab at the Georgia Institute of Technology.

1.        R. Ponzini et al., IEEE Trans. Biomed. Eng. 57, 1807–1815 (2010).

2.        I. C. Campbell et al., J. Biomech. Eng. 134, 51001 (2012).

Introduction
An obstacle in modeling is the estimation of consistent parameters, for which outputs of the simulation reproduce some given data. In particular when the number of parameters exceed the number of outputs, all parameters cannot be estimated uniquely, increasing the uncertainty. We tackle this problem here through MAP estimation, which computes the most likely parameters according to a prior distribution. We compare two types of estimation on a complex parameter estimation problem of 13 parameters of a 0D cardiovascular model from data of 11 patients with Pulmonary Hypertension before and after therapy.

Methods
We build a 0D model of the cardiovascular circulation, based on the 0D model of ventricular mechanics introduced in [1]. This cardiovascular model has 4 connected components: the left heart, the right heart, the systemic circulation and the pulmonary circulation. It is implemented in the language CellML and will be made available, easily exploited through the software OpenCOR.

We then estimated, for 11 patients and at two instants (22 cases), 13 cardiovascular parameters (resistances, compliances, contractilities/radiuses of the ventricles) from 10 outputs (pressures values and volume), through the optimisation of a regularised cost function with a genetic algorithm. We compare a parameter estimation without priors, and a MAP estimation with "population-based" priors. Namely, we use as the prior distribution for this second estimation the distribution of the population of estimated parameters without priors.

Results
First, in both parameter estimation, all the 22 but one simulations have their 10 outputs fit almost perfectly the target data, demonstrating the efficiency of the optimisation framework on such complex problem. We also observe that the use of priors in MAP estimation promotes small variation of all parameters and limits large variations of few parameters when multiple parameters values are possible, which is more consistent for biophysical values.

Secondly, we show that the estimated parameters of resistance in the model are very correlated to the traditionally estimated resistance in clinical practice (from the pressure values), demonstrating consistency of the model with respect to clinical practice. In particular, variation of estimated parameters before and after therapy reflects clearly which patients are non-responders to the therapy.

Discussion
We present an efficient framework to build personalised simulations of the cardiovascular system and compare their parameters. Beside the specific application to Pulmonary Hypertension we believe out framework is flexible and could be used for many other problems. The population-based priors enable to constrain the estimation in order to reduce the uncertainty coming from the limited observability of the parameters based on clinical data.

References
[1] Molléro et al. Multifidelity-CMA: A Multifidelity Approach for Efficient Personalisation of 3D Cardiac Electromechanical Models. Biomechanics and Modeling in Mechanobiology, Sep 2017.

Introduction
For decades fluid dynamics models have been used to predict flow and pressure wave propagation in the cardiovascular system [1]. Recently focus has shifted to fitting measured waveforms using geometries from imaging [2]. Yet little work has been done to understand how error in image segmentation impacts fluid dynamics. This is of particular importance in network models, since at each bifurcation the vessel dimensions change significantly. In mice pulmonary vasculature, radii change rapidly with generation, e.g. the main pulmonary artery has radius > 0.5 mm while vessels downstream have significantly smaller radii. Additionally, network connectivity varies from specimen to specimen. This study addresses uncertainty in blood flow and pressure predictions using a 1D network extracted from micro-CT (μCT) images together with flow and pressure measurements from control mice.

Methods
1D Networks are extracted from μCT images from C57BL6/J mice (Jackson Laboratory, Bar Harbor, ME) using global thresholding techniques. Lower threshold and smoothing image segmentation parameters are varied for a single image. For our analysis, we fix the number of steps for a semi-automated contour evolution algorithm. The 1D networks are a collection of centerlines connected at the barycenter of each junction. Typical networks analyzed contain 40 or more vessels.
All vessel segments are assumed straight with a radius defined as the mean over the vessel, calculated away from the junction. Fluid dynamic predictions are obtained by solving the 1D Navier-Stokes equations combined with a linear constitutive equation relating pressure and area (for details see [2-4]). At junctions, pressure is assumed continuous and flow is conserved. A measured flow profile is applied at the inlet and a three element Windkessel boundary condition is attached at each terminal vessel. Similar to previous studies [3] nominal parameter values were calculated using imaging data and available flow and pressure data.

Results and Discussion
Figure 1 shows an example network along with pressure and flow computations from 1000 simulations. Each simulation is computed using a variable geometry, where the variance of each vessel’s dimensions comes from 20 a priori segmentations of one μCT image. A measured inflow is attached at the inlet, so the main variations are associated with predictions of pressure and vessel area (not shown). Results reveal significant variation within the network, emphasizing model sensitivity to geometric parameters. Uncertainty is also associated with other model parameters including with outflow boundary conditions, vessel stiffness, and the inflow profile [4].

Acknowledgement
NSF-DMS-1615829

 References
[1] van de Vosse et al. (2011) Annu Rev Fluid Mech, 43:p467
[2] Qureshi et al, (2014) Biomech Model Mechanobiol 13: p1137
[3] Qureshi et al. (2017) arXiv:1712.01699 [physics.flu-dyn]
[4] Paun et al. (2017) Stat Neerl, in press.

Introduction
For clinical relevance of simulation of cardiovascular processes and diseases, it is essential to provide more than just the result of one deterministic simulation, even if this is cutting edge and is incorporating all sorts of complexities. One also needs to provide for example information about parameter sensitivity of results and sort of quantify the overall level of confidence. This means that instead of using population mean values to perform patient-specific simulations, neglecting inter- and intra-patient variations present in these parameters, such uncertainties have to be considered in the computations. This can be achieved by including uncertainties of input quantities directly in the simulation and by calculating the stochastic information of the quantities of interest. However, due to limited computational resources and several shortcomings of traditional UQ approaches, parametric uncertainties, modeled as random fields, have not yet been considered in patient-specific, nonlinear, large-scale, and complex biomechanical applications.

Methods
We present an efficient UQ scheme able to deal with complex large-scale problems in biomechanics. The UQ framework is based on multi-fidelity sampling and Bayesian formulations. The key feature of the presented method is the ability to rigorously exploit and incorporate information from low fidelity models. Most importantly, response statistics of the corresponding high fidelity model can be computed accurately even if the low fidelity model provides only a very poor approximation of the quantities of interest. The approach merely requires that the low fidelity model and the corresponding high fidelity model share a similar stochastic structure. This results in a tremendous flexibility in choice of the approximate models.

Results & Discussion
Flexibility and capabilities are demonstrated via UQ for some challenging problems. We start with patient-specific, large-scale, nonlinear FE models of AAAs. Here, constitutive parameter of the wall are modeled as univariate 3D, non-Gaussian random fields, thereby taking into account inter-patient as well as intra-patient variations. Direct Monte Carlo demonstrates excellent quality of results. Additionally, the employed approach results in a tremendous reduction of computational costs, rendering UQ with complex patient-specific nonlinear biomechanical models practical for the first time. We will then show the general applicability of the approach via applying it to UQ for cardiac simulations, using low fidelity models obtained from reduced order modeling. If time allows we will also apply the approach to respiratory biomechanics.

References
Biehler J., Gee M.W., Wall W.A.: Towards efficient uncertainty quantification in complex and large scale biomechanical problems based on a Bayesian multi fidelity scheme. Biomechanics and Modeling in Mechanobiology, 14(3) (2015), 489-513
Biehler J., Kehl S., Gee M.W., Tanios F., Pelisek J., Maier A., Reeps C., Eckstein H.-H., Wall W.A.: Non-Invasive Prediction of Wall Properties of Abdominal Aortic Aneurysms Using Bayesian Regression. Biomechanics and Modeling in Mechanobiology, 16 (2017), 45-61

We present new trends for parametric computing based on advanced offline-online POD-Galerkin reduced order methods for applications in computational fluid dynamics, with a special interest in problems related with optimisation, control, uncertainty quantification, as well as data assimilation, for parametric cardiovascular flows. These problems are of particular interest due to the fact that they bring a certain complexity related with the geometry of biological structures to be reconstructed and parametrised for the simulation, as well as patient specific medical data subject to uncertainty and noise. In addition the numerical approximation of complex nonlinear parametric systems needs proper stabilisation also at the reduced order level. Few examples from real-life problems are introduced as proof-of-concept to present a complete computational pipeline in an advanced reduced order setting. Multi-physics is accounted in considering fluid and structure interaction between flows and arterial walls, and optimal flow control is accounted for the solution of inverse problems (defective boundary conditions, parameter estimation, data assimilation) as well as for interfacing network components. A related topic of interest is the detection of flow bifurcations in cardiovascular settings, as well as Coanda effect. This work is in collaboration with Dr Laura Jimenez at Sunnybrook Hospital in Toronto (Canada), Prof. Piero Triverio at University of Toronto (Canada) and Prof. Annalisa Quaini at University of Houston (USA).

Introduction
Hemodynamic forces have been shown to affect progression of cerebral aneurysms.¹ Blood flow dynamics can be assessed in vivo, using 4D Flow MRI or with patient-specific computational fluid dynamics (CFD). The dynamic range of the 4D Flow is determined by a velocity sensitivity parameter (venc), set above the expected maximum velocity, which can result in noisy data for slow flow regions. A dual-venc 4D Flow MRI technique, where both low- and high-venc data are acquired, can improve velocity-to-noise ratio and, therefore, quantification of clinically-relevant hemodynamic metric.² Alternatively, CFD provides superior spatiotemporal resolution, but its accuracy depends on modelling assumptions. In this study, we compare single- and dual-venc in vivo 4D Flow MRI to CFD results to determine the advantages of dual-venc relative to a single-venc acquisition.

Methods
4D Flow and MR angiography (MRA) data were acquired for four cerebral aneurysm patients at Northwestern University School of Medicine. The MRA datasets were segmented to generate patient-specific geometries. Navier-Stokes equations were solved with the finite-volume solver Fluent (ANSYS). The inlet and outlet waveforms were prescribed from the 4D Flow measurements. The flow fields obtained with CFD and both single- and dual-venc were compared to determine differences between techniques and modalities. A cross-sectional plane was placed within the aneurysm geometry to show both high-velocity jets and low-velocity vortices.

Results
One case is presented here due to space limitations. The 4D Flow and CFD results obtained for an anterior communicating artery (Acom) aneurysm are shown in Figure 1. Panels (A)-(C) display flow streamlines obtained with the single-venc and dual-venc 4D Flow, and with CFD. The dual-venc acquisition enables visualization of the slow, recirculating flow in the aneurysm. Panels (D) and (E) display the difference between single-venc and dual-venc with CFD data, respectively. Panel (F) compares single- to dual-venc measurements.


Figure 1: Flow streamlines obtained with A) single-venc; B) dual-venc 4D Flow; C) CFD; and velocity difference maps between D) single-venc and CFD, E) dual-venc and CFD, and F) single- and dual-venc 4D Flow.

Discussion
The dual-venc technique provides a greater range of velocities; however, MRI velocimetry still may have inadequate resolution for estimation of clinically-relevant metric. Velocities acquired with dual-venc were 2.51% better matched to CFD within the aneurysm. The difference maps show only slight deviation between the single- and dual-venc velocities when compared to CFD with averaged absolute error giving: (D) 0.1049, (E) 0.0932, and (F) 0.0436. The preliminary data in other patients indicate that dual-venc 4D Flow MRI provides similar enhancement of velocities; further quantification of the errors on all available data is necessary.

References
1. Schnell, S., et al., (2017). J Magn Reson Imaging, 46(1) p102

Introduction

The QT interval of an electrocardiogram signal, which begins with the activation of the ventricles and ends with their recovery, is a function of both heart rate and ventricular repolarization time. A prolonged QT interval is typically indicative of a patient’s vulnerability to Torsade de Pointes, a potentially fatal arrhythmia characterized by rapid depolarization and repolarization of the ventricles. Various drugs have undesired side effects associated with a prolongation of the QT interval and an increased pro-arrhythmic risk. Two landmark parameters that quantify the effect of the drug are the drug concentration IC50 at 50% blockage of an ion channel, and the slope, h of the dose-response curve at this point. 

Methods

Here we present a novel high resolution, multi-scale, multi-fidelity computational model to predict the QT interval for a wide range of IC50 and h values [1]. The experimental variability of these two parameters is the origin of uncertainty in the model. Leveraging recent developments in hierarchical Bayesian inference [2], we sample a set of 500 IC50 and h parameter pairs from the experimental dose-response data for 30 common drugs [3]. For each pair of parameters and a chosen drug concentration, we calculate the conductance block of a specific ion channel. Using these ion channel blocks, we construct a probabilistic multi-fidelity model to predict a set of QT intervals using Gaussian process regression. Finally, we employ Gaussian kernel density estimations to produce the probability density function of the QT interval at each concentration.

Results

Our uncertainty quantification reveals that, for a compound with three dose-response experiments, the QT interval has up to 76% variability at the maximum concentration compared to a baseline case of zero concentration and zero ion channel blockage.

Discussion

We have quantified the effect of experimental errors in dose-response curves in a critical biomarker of drug-induced arrhythmias: the QT interval. Our results showcase a novel and efficient workflow for predicting the pro-arrhythmic risk of common drugs and assessing their effect on cardiovascular performance.

References

1. F. Sahli, J. Yao, E. Kuhl. Predicting drug-induced arrhythmias by multiscale modeling, submitted (2017)
2. R. H. Johnstone, R. Bardenet, D. J. Gavaghan, G. R. Mirams, Hierarchical bayesian inference for ion channel screening dose-response data, Wellcome open research 1, 6 (2016).
3. W. J. Crumb, J. Vicente, L. Johannesen, D. G. Strauss, An evaluation of 30 clinical drugs against the comprehensive in vitro proarrhythmia assay (cipa) proposed ion channel panel, Journal of pharmacological and toxicological methods 81 (2016) 251–262. 


Computational fluid dynamics models are increasingly used to simulate vascular flow. Patient-specific flow measurements are typically used to set flow boundary conditions that drive the simulations. When patient-specific flow measurements are unavailable, mean values of flow measurements across small cohorts are used as normative values. In reality, both the between-subjects and within-subject flow variability is large, and mean values across a cohort might not be indicative of flow in a particular individual. We highlight the importance of quantifying the uncertainty in CFD-based indicators, such as wall shear stress (WSS), due to the physiological variability of blood flow.

We develop data-driven statistical models for the between-subjects variability of internal carotid flow (Lassila et al. 2018). A log-linear mixed effects model describes the variability in the time-averaged flow rate, while a Gaussian process model describes the variability of the time-dependent flow fluctuations. The model parameters are identified from carotid ultrasound measurements in a cohort of 104 elderly volunteers. To model within-subject variability, we use a cerebral autoregulation model (Mader et al. 2015) to model the response of internal carotid flow to changes in heart rate and blood pressure. Our models provide: (i) cohort-specific boundary conditions to vascular CFD simulations in the case that patient-specific flow measurements are not available, and (ii) estimates of within-subject flow variability in cases where patient-specific flow measurements are available only as spot measurements.

The flow variability models are applied to study the effects of flow variability on WSS patterns in intracranial aneurysms. Different WSS-related indicators have been proposed for identifying aneurysms with high rupture risk. We have previously shown that two commonly used WSS-based indicators, oscillatory shear index (OSI) and transverse WSS (TransWSS), are sensitive to physiological flow variability (Sarrami-Foroushani et al. 2016), which may confound their interpretation. We further study the effects of exercise-induced variability of ICA flow on WSS-based indicators. While the cerebral autoregulation system holds the time-averaged flow rate relatively stable even under strenuous exercise, changes in the flow waveform are more dramatic and can lead to large uncertainty in WSS-based indicators.

A thoracic aortic aneurysm is a typical cardiovascular disease consisting in a permanent dilatation of the thoracic aorta, and thus represents a symptom of weakness of the arterial wall. Even though the risk of serious complications depends on the complex biomechanical state of the artery, clinical strategies are still based on the simplistic criterion of the maximum aortic diameter, which has been shown to be far from adequate, especially for small and medium-sized aneurysms.

Since the literature has pointed out the decisive impact of hemodynamic indicators in the origin and evolution of aneurysms, computational fluid dynamics has been extensively used in the last few years. Numerical simulations allow a large number of variables that cannot be measured in-vivo with sufficient precision to be obtained. However, the accuracy of the results strongly depends on the assumptions made to build the hemodynamic model.

We focus herein on inlet boundary conditions based on the imposition of flow rate waveforms. We aim to assess the impact that uncertainties in inlet data have on the hemodynamic descriptors. This is motivated by the fact that patient-specific measured data are not always available. We also wish to highlight how possible variations of the different quantities characterizing the inflow waveform affect the hemodynamic indicators.

A patient-specific geometry of thoracic aortic aneurysm was acquired by means of Magnetic Resonance Imaging (MRI), together with the relevant inlet flow rate waveform, which was taken as a reference for the subsequent stochastic analysis. The selected uncertain parameters were the stroke volume and the cardiac cycle period. We first assumed uniform PDF distributions, with specified variation intervals of the parameters, then more accurate beta PDF distributions obtained from clinical data about 23 patients were used. All the analyses were performed with the generalized Polynomial Chaos approach and by using the open-source code SimVascular for the single deterministic simulations, which were carried out for both rigid and deformable walls.

In general, we found a large influence of the two input parameters on Time-Averaged Wall Shear Stress (TAWSS). In the case of clinically-based beta PDFs, the space-averaged TAWSS stochastic standard deviation is about 0.4 Pa (almost 50% of the space-averaged reference deterministic value) and the stroke volume shows a clear predominant effect all over the arterial wall. On the other hand, the two parameters have a comparable impact on TAWSS in case of uniform-assumed PDFs. Finally, the wall compliance seems to have a limited overall effect: indeed, in the deterministic simulations, the wall shear stress distributions predicted by the rigid and deformable models are very similar, especially for large values of stroke volume.

Introduction
The parameters that characterize reduced 1D models of the cardiovascular system (CVS) feature variability both between patients and within a single individual. Characterizing this variability via uncertainty quantification (UQ) studies on the CVS entails two major challenges: (i) The computational resources required by full 3D models may not be easily accessible by users (e.g., hospitals); (ii) Reduced 1D models may be inaccurate in capturing anomalies of the physiology in the presence of cardiovascular pathologies. This introduces additional (epistemic) uncertainty that should also be characterized.
The objectives of this study are (I) to design UQ solvers tailored to the CVS that promote parallelism and scalability, and (II) to enhance the accuracy of 1D reduced models without incurring computational cost.

Methods
The Domain Decomposition Uncertainty Quantification (DDUQ) approach [1] performs UQ at the subsystem level, and propagates uncertainty information encoded as polynomial chaos coefficients via overlapping DD techniques, allowing for a reduction of the computational time.
The local cross-sectional dynamics discarded by 1D models can be retrieved via educated reduced models such as the Transversally Enriched Pipe Element Method (TEPEM) [2], or HiMod methods [3]. The main stream dynamics and the transverse components are solved by the Finite Element Method (FEM) and by Spectral Methods, respectively, to guarantee high accuracy at low computational cost.

Results
Preliminary results obtained for a 2D steady non-linear heat equation highlight the scalability of the DDUQ approach, with different smoothers and relaxation parameters (Fig. 1 - left) [1].
The convergence rate of the TEPEM solver for the solution of the unsteady Navier-Stokes equations in a 3D tortuous pipe is by far higher than for FEM (Fig. 1 - right), with a reduction of the computational time up to 900 times [2].


Discussion
This study shows that (i) The computational cost of UQ in the CVS can be drastically reduced by promoting the independence of the subsystems and avoiding full-system simulations; (ii) The solver is scalable; (iii) Educated reduced models improve the accuracy of 1D solvers, at approximately the same computational cost.

Acknowledgements
CNPq, FAPERJ, NSF grant DMS 1419060.
Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525.

References

1. K. Carlberg, M. Khalil, K. Sargsyan, and S. Guzzetti. Uncertainty propagation in large-scale networks via Domain Decomposition. In preparation, 2017.
2. L. Mansilla Alvarez, P. Blanco, C. Bulant, E. Dari, A. Veneziani, and R. Feijóo. Transversally enriched pipe element method (TEPEM): An effective numerical approach for blood flow modeling. IJNMBE, 33(4), 2017.
3. S. Guzzetti, S. Perotto, A. Veneziani. Hierarchical Model Reduction for Incompressible Fluids in Pipes. IJNME. 2017.