Physics driven reduced order model for real time blood flow simulations

11/05/2019
by   Sethuraman Sankaran, et al.
0

Predictive modeling of blood flow and pressure have numerous applications ranging from non-invasive assessment of functional significance of disease to planning invasive procedures. While several such predictive modeling techniques have been proposed, their use in the clinic has been limited due in part to the significant time required to perform virtual interventions and compute the resultant changes in hemodynamic conditions. We propose a fast hemodynamic assessment method based on first constructing an exploration space of geometries, tailored to each patient, and subsequently building a physics driven reduced order model in this space. We demonstrate that this method can predict fractional flow reserve derived from coronary computed tomography angiography in response to changes to a patient-specific lumen geometry in real time while achieving high accuracy when compared to computational fluid dynamics simulations. We validated this method on over 1300 patients that received a coronary CT scan and demonstrated a correlation coefficient of 0.98 with an error of 0.005 +- 0.015 (95 compared to three-dimensional blood flow calculations.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 6

page 10

page 13

page 14

11/08/2021

Automated generation of 0D and 1D reduced-order models of patient-specific blood flow

Three-dimensional (3D) cardiovascular fluid dynamics simulations typical...
05/29/2018

Comparison of 1D and 3D Models for the Estimation of Fractional Flow Reserve

In this work we propose to validate the predictive capabilities of one-d...
07/06/2020

Non-intrusive PODI-ROM for patient-specific aortic blood flow in presence of a LVAD device

Left ventricular assist devices (LVADs) are used to provide haemodynamic...
03/29/2021

Visualizing Carotid Blood Flow Simulations for Stroke Prevention

In this work, we investigate how concepts from medical flow visualizatio...
12/04/2021

Modeling and Predicting Blood Flow Characteristics through Double Stenosed Artery from CFD simulation using Deep Learning Models

Establishing patient-specific finite element analysis (FEA) models for c...
01/24/2019

Periodic-corrected data driven coupling of blood flow and vessel wall for virtual surgery

Fast and realistic coupling of blood flow and vessel wall is of great im...
01/05/2022

Nonlinear lumped-parameter models for blood flow simulations in networks of vessels

To address the issue of computational efficiency related to the modellin...
This week in AI

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

1 Introduction

Patient-specific modeling of blood flow has emerged as a tool with increasing importance in the diagnosis and treatment of patients with coronary artery disease [37] [28] [3] [13] [24] [14]. A unique strength of simulation methodologies lie in the predictive modeling of hemodynamics in response to unplanned events (such as progression or regression of lesions), or the outcome of planned procedures (such as surgical intervention) [36]. For predictive modeling to be used in the cardiac catheterization laboratory or other invasive procedures, it is imperative that the modeling tools can generate results in seconds.

In this work, we address one of the most challenging aspects of interventional planning for patients with coronary artery disease - Which coronary artery stenoses are having the greatest impact on blood flow? Due to the vastly different anatomies and lesion morphologies in a normal patient population, it is infeasible to have a preset approach for the treatment strategy. Therefore, the planning of invasive procedures in a patient is left to the knowledge, intuition and experience of the clinician based on available data. This paper describes a patient-specific framework that can virtually model different scenarios and help the physician evaluate various strategies by identifying and ranking the impact each coronary artery stenosis has on the blood flow. We describe a fast and accurate tool derived from a patient specific coronary CT angiography scan that is applicable in a clinical setting. The clinical metric computed is the fractional flow reserve (FFR) which is measured in the cardiac catheterization lab with a pressure wire inserted in the patient’s coronary arteries. FFR is ratio of the time-averaged pressure downstream of a coronary artery stenosis to a time-averaged reference aortic pressure under conditions of maximum hyperemia typically induced by the intravenous administration of adenosine. FFR computed from coronary computed tomography data () is derived by simulating blood flow and pressure in a patient-specific anatomic model using computational fluid dynamics (CFD) methods [37].

Fast assessment of blood flow has been an active area of research in the past several years. Recently, many successful methods have evolved that apply simulation methodologies for the estimation of clinical quantities of interest from medical imaging data. For instance, artificial intelligence methods have been applied to quantify hemodynamics from phase contrast MRI scans 

[25], or to quantify information from handheld ultrasound devices [5]. While these techniques help in the automation, miniaturization or reproducibility in the extraction of information, the approach, fundamentally, is to quantify information already present in the medical images. In contrast, this paper focuses on predictive modeling of the effect of treatment, for which no current patient-specific non-invasive alternatives exist.

The development of reduced order models for the simulation of blood flow has been the focus of several studies [34, 4, 19, 27, 18, 11, 15]. Kassab and colleagues [20] proposed an analytical model based on the conservation of energy, and considered various energy losses associated with a lumen narrowing such as convection and diffusion and also accounted for the losses associated with sudden constriction and expansion in lumen area. They demonstrated the performance on a tube-like stenosis model, however calculation of some of the higher order geometric quantities (higher order gradients of radius or area) accurately on patient-specific geometries is not possible due to the imaging resolution. Schrauwen et al. [9] [10] improved upon the traditional reduced order models that assume parabolic velocity by deriving a velocity profile based on the geometry and flow, and validated it on straightened coronary arteries. However, the velocity profiles on patient specific models with multiple bifurcations and interfering stenosis (effect of one stenosis overlapping on another) can be significantly different. Nithiarasu and colleagues [1] [7] [8] developed reduced order models that accounted for intramyocardial pressure and material properties of the arterial wall, and evaluated performance in a virtual cohort of 30 lesions with the same global geometry. Itu et al. [12]

proposed a machine learning approach for assessment of fractional flow reserve and verified it against another reduced order model on synthetic data. While these approaches have helped immensely in understanding how the complex physics of fluids can be captured with a reduced order system, these are not applicable for clinical translation due to the high demands on accuracy, latency and robustness. Our goal in this paper is to develop an algorithm that meets these requirements by demonstrating a very high accuracy over a large cohort of patient-specific models with diverse charecteristics.

The framework we propose uses a response surface [22] methodology and can predict the results of the simulations significantly faster and with an accuracy close to the high-fidelity simulations. The reduced order model is parameterized using the response surface built using a full order high fidelity high cost (HFHC) model performed at certain pre-set configurations. Simulations for the HFHC model may utilize all the information available about the system (such as using the full spatial and temporal representation). A response surface is a mathematical relationship between a quantity or quantities of interest or parameters, and an underlying representation or property, and have been built successfully for optimization problems [32]

. The pre-set configurations involve exploration of both the space of lumen geometries and flow rates. As we increase the number of HFHC simulations used to construct the response surface, the accuracy of the real-time predictive model becomes closer to the HFHC model within some error margin. The response surface is used to explore the parameteric space by using a reduced order model to interpolate the solutions. As a result of using physics driven interpolatory functions derived from one-dimensional Navier-Stokes equations, our approach can mimic the results of the computational fluid dynamics simulations better than other standard approaches such as using polynomial or Lagrange interpolatory functions. The response surface is built to exactly match the output of the high-fidelity model for the configurations where high-fidelity simulations were performed.

We demonstrate that we can achieve excellent accuracy of the proposed method compared to three-dimensional simulations for coronary blood flow using just four full order offline simulations. We also evaluated the performance of the algorithm in predicting hemodynamics post stenting or percutaneous coronary intervention (PCI) against invasive data for different lesion types. PCI is the process of inserting a catheter with a balloon attached to restore and attain ideal lumen geometry on a section of the blood vessel. We demonstrate that we can achieve high accuracy against invasive data and that it can be used as an interactive tool to plan PCI.

We structure the paper as follows. In section , we describe the method behind the patient-specific reduced order model. We describe the HFHC model and demonstrate how we can build a real-time model using a set of HFHC models. In section , we demonstrate the performance of the method against HFHC simulations. We use internal data that encompasses different lesion types and FFR ranges to validate the model. We finally discuss implications of the method and clinical translation in section .

2 Methods

Blood flow in the human arterial system can be effectively modeled using computational fluid dynamics [37] [35]. These simulations can be performed on complex patient-specific anatomic geometries that are reconstructed from coronary CT scans. The Navier-Stokes equations are used to model blood flow and lumped parameter boundary conditions are prescribed that account for the flow-pressure relationship in the microvasculature (vessels that are not visible in the CCTA and therefore not modeled). Vignon-Clementel et al. [16] developed methods for coupling boundary conditions using a circuit analogy of the microvasculature with the Navier-Stokes equations and Kim et al. [23]

specialized this method to modeling boundary conditions for coronary arteries. Sankaran et al. used a set of ordinary differential equations 

[17] instead of an analytical formulation for modeling blood flow in the microvasculature. In this work, we use resistance (ratio of pressure drop to flow) boundary conditions for the microvascular network, which are now routinely used in clinical practice [37] for the non-invasive assessment of .

2.1 Computational Fluid Dynamics Model

The partial differential equations that model blood flow are represented as


(1)

with the boundary conditions

(2)

where is a non-linear differential operator modeling the Navier-Stokes equations, are the blood velocities and pressures, represents fixed parameters, is the problem domain and is the boundary of the domain. Since this model is already validated [28] [14] and used in clinical practice for the calculation of , we use the model governed by Equations 1 and 2 as our HFHC model. Our goal in this work is to build a low cost model that is equivalent in accuracy to the high fidelity model.

2.2 Reduced Order Models

A reduced order model of the partial differential equation approximates the operator, , using a simpler operator, (e.g., ordinary differential equations), reduces the dimensionality of the solution space and/or reduces the dimensionality, of the problem space , to a lower dimensional space, , and simplifies the parameter set to . In this case, can be the three dimensional representation of the geometry, and can be a corresponding one-dimensional representation. The reduced order model problem is posed as follows, where we omit the spatial dependence of for clarity:

(3)

with the boundary conditions

(4)

The goal is to have be a good approximation of at . To achieve this, our approach involves performing simulations of the original HFHC model for various geometries, , boundary domains, and boundary conditions such that we can generate a response surface and subsequently build a good approximating solution to the problem. In order to generate the response surface and to ensure that we limit the exploration space to realistic geometries, we can impose bounds on the domain of interest and the boundary conditions that the system will be subject to. The original governing equation is solved in a series of domains, boundary conditions and parameters.

(5)

where is the number of HFHC simulations performed. The results of the HFHC simulation for each of these configurations are

(6)
Figure 1: (left) The reduced order model has two terms: a purely geometric term (parameterized by radius ) and a term that varies linearly with flowrate, , with a geometry dependent slope, . The coefficients are fit for each patient based on the patient geometry and idealized geometry, and (right) four simulations are performed on two lumen segmentations (patient and idealized) with two boundary conditions provided for each.

2.3 Physics driven response surface method

Any sampling method or quadrature method can be used for the selection of the configurations in Equation 5. As each simulation is computationally expensive, sampling methods such as Monte-Carlo or latin hypercube sampling [26] tend to converge slowly in building an accurate representation of the underlying function. Functional space methods, such as stochastic collocation or adaptive stochastic collocation methods [38] [33], can help accelerate the convergence for a class of functions. However, they still demand significant calculation to generate the response surface. We improve upon this further by regularizing the problem using a physics driven interpolation approach. Instead of using polynomials, which is commonly used for response surfaces, we use the analytical solution from one-dimensional, steady state Navier-Stokes equations to interpolate within the response surface. This enables employing a significantly reduced number of simulations while still enabling equivalence with the HFHC simulations. In fact, we later demonstrate that four HFHC simulations at the extremas of the exploration space are sufficient to obtain equivalence with the full order model (refer Fig. 1).

A patient-specific response surface, , is a mapping of the coefficients of the resistance given a geometry, boundary conditions and parameters,

(7)

wherein captures the complexity of the original equations, enabling to be a less complex operator than . The reduced order model, , can be constructed such that at the configurations where the HFHC simulations are performed. This approach lets us solve the reduced order problem (refer to equations 3 and 4) for the unknowns , while ensuring that the results are identical at the configurations. The solution field contains only flowrates and pressures instead of three velocity components and pressures. The reduced model is evaluated only along the vessel centerlines along the vessel path instead of the entire volumetric patient-specific geometry. The reduced order operator contains the affine resistance model that relates flowrates to pressures and the flow-split model at bifurcations, which are explained in later sections. Naturally, the approximations to the HFHC solutions at the intermediate configurations will be better for larger , but so will the time needed for the offline computations. These are described in detail in the following sections.

Figure 2: Illustration of the proposed method that builds a patient-specific reduced order model which depends on both the anatomy and the physiologic state. Four simulations are used to explore the embedding space of the lumen geometry and physiologic state. Hyperemia and Superemia refer to different boundary conditions (or flowrates) and the colormap within the window represents how the reduced order model resistance varies within the patient-specific exploration space. is the local pressure and is the reference aortic pressure.

2.4 Patient-specific ideal geometries

One of the novelties of our work is to build a patient-specific exploration space. To define the bounds on the geometric domain, at each location in the coronary tree, we consider modifications wherein the vessel radius falls between the patient radius and the idealized radius. The idealized model represents a model of a patient without regions of lumen narrowing (or “healthy looking”). The idealized model is constructed by seeking a radius profile that is monotonically non-increasing from the ostium to the leaves and is closest to the original radius profile, using a robust objective function, formulated mathematically as

(8)

such that

(9)
(10)

where is the input radius, is the idealized radius estimate and is an optional minimal radius constraint at a vessel tree location i. For instance, can be set equal to to force the ideal model to be point-wise at least as large as the patient model. Idealized radius values are optimized globally, over the entire coronary tree, through a dynamic programming approach. Given idealized radius values, a patient-specific idealized 3D model is constructed by dilating the original model from its local radius to the target idealized radius . This step is performed by merging 3D spheres of radii into the implicit representation of the original geometry (signed distance field) and generating a new mesh from that dilated implicit representation. Figure 3 shows an ideal model (and radius profile) superimposed over the patient model. We omit the index to represent the entire lumen segmentation of the appropriate radii.

Figure 3: (left) A section of a vessel with patient lumen segmentation (solid) and idealized lumen segmentation (meshed) and (right) the corresponding patient radius (blue) and idealized patient radius (dotted orange).

2.5 Flow dependent hemodynamic resistance

Lumen geometry influences the flowrate through the model. Therefore, quantifying the impact of geometry on hemodynamics should account for the changes in the flowrate through the model. The patient-specific reduced order model (PSROM) assumes an affine resistance model for flow through the coronary arteries,

(11)

where is the local resistance of the lumen segmentation to blood flow, is the blood flow rate, is the resistance to flow that is purely geometric, and is the sensitivity of resistance to flowrate. This is a generalization of having a viscous (Poiseuille) and inertial (Bernoulli) term, but the coefficients being derived from HFHC simulations instead of assuming a circular tube cross-section and a parabolic velocity profile [30] [31].

To accurately build a response surface for and , it is imperative that we explore the space of flowrates in the patient-specific model. Since flowrate is not an independent variable, but rather depends on the boundary conditions and the vessel lumen geometry, we modify the boundary conditions to change the state of the flowrate to derive . We parameterize the lumen geometry by its radius, which is commonly used for analytical 1-dimensional modeling of the Navier-Stokes equations [30]. This modification can be done for both the patient-specific geometry and the idealized lumen geometry, and interpolated for intermediate configurations.

Since we have chosen to use an affine model, two flowrates should suffice to explore the space of flowrates. Since the patient-specific model is already solved for the hyperemic state, we need only one additional flowrate. In this paper, we are interested only in dilation of the lumen geometry (not constriction) which would result only in the increase in the rate of blood flow. Therefore, we consider an additional state, called henceforth as “superemia”, wherein we reduce the microvascular resistance by 40%. This number was chosen based on a grid search on a development dataset. Smaller values such as 10% or 20% tend to have less signal to noise (low flow difference and pressure loss compared to the numerical noise of the CFD simulations). Larger values such as 80% tend to have flowrates very different from the modified lumen geometries. An illustration of the hyperemia and superemia configurations on patient and idealized geometrieis can be seen in Figure 2.

2.6 Patient-specific Reduced order model

The patient-specific reduced order model (PSROM) is constructed by first calculating and for the original and idealized models, where two simulations are performed for each model (hyperemia and superemia). For any new configuration, the intercept can be calculated as

(12)

where

(13)
(14)

and

(15)

In addition, the slope can be calculated as

(16)
(17)
(18)

and

(19)

and are the patient and idealized lumen areas, and derived directly from the radii. The interpolation functions, and were chosen to be consistent with the Poiseuille pressure loss and Bernoulli pressure loss respectively. A predictor corrector algorithm [2] is used to solve the equations governing PSROM. Further details of how these quantities are calculated are given in Algorithm 1. Figure 4 shows a schematic of the proposed patient-specific ROM approach.

2.7 Predictor corrector algorithm

The goal of the proposed method is to predict hemodynamics given a geometry modified from the original patient’s geometry. For the modified geometry, the response surface that is constructed from the HFHC simulations is used to estimate the coefficients and given a radius . These are then used to calculate the local resistances given a flowrate . These resistances are integrated from the leaves of the coronary tree using a circuit analogy (resistances are solved in series within a given segment, and in parallel at each bifurcation). The resistances at the ostium are then used to calculate an updated flowrate. These steps are repeated until a pre-set convergence criterion is met. The boundary conditions depend on the hemodynamics, and we account for this dependence by modeling the average size of the downstream vasculature being proportional to the pressure at the outlet (passive response to increased pressure). At each outlet, the boundary conditions are scaled from the original value by a ratio of the outlet pressures.

There are two conditions where we do not apply the interpolation approach. In regions just downstream of the region of lumen modification, the original patient-specific geometry may have negative resistance (increase in pressure) and thereby a negative slope. Such regions are absent when the initial section with lumen narrowing is replaced by an idealized section. We empirically define this potential pressure recovery region to be up to distal to a region of vessel narrowing. This was also chosen based on the development dataset. We might miss capturing pressure increase if we choose too small a pressure recovery region length cutoff. On the other hand, the choice of larger region lengths may result in inclusion of recovery that are not related to the flow features associated with a stenosis. Further, we assume that pressure can recover only if the area gradient is negative. We therefore replace any region with negative area gradient in the pressure recovery region with a Poiseulle resistance model. In addition, we do not use the slope and intercept estimated from the four solutions if the hyperemic and superemic flowrates are too close to each other. This is because numerical convergence related errors in the HFHC model can have a non-linear effect on the slope and the intercept, and the signal to noise effect can result in lower accuracy of the prediction, especially when flowrates of the modified geometry are significantly larger than the original patient-specific model.

Figure 4: (A) Steps involved in generating the patient specific reduced order model. First, a lumen segmentation of the patient anatomy is constructed from a coronary CT scan. Next, an idealized model of the patient’s coronary tree is constructed that represents a segmentation without lumen narrowing. Two sections of the model are highlighted to illustrate differences between the patient and idealized models. Next, four CFD simulations are performed, two each on the patient and idealized models with two different boundary conditions. (B) the patient specific reduced order model (PSROM) can be used on any modified geometry that lies between the patient and idealized geometry to calculate the hemodynamics.

2.8 Overall Algorithm

The overall algorithm is given in Table 1 and illustrated in Figure 4. A convergence criteria of for the ostial flowrate is used. The extraction of radii from the patient-specific lumen segmentation and the subsequent generation of the idealized model are performed prior to building the physics driven response surface. represents the set of all centerline points, represents the set of forward neighbors (implying leaf centerline points do not have a neighbor), is the pressure at centerline point , is the flowrate of the section connecting i to k, and are the slope and intercept of the section connecting centerline to , superscripts and represent superemia and hyperemia respectively. For the predictor corrector method, the superscripts denote the iteration number. denotes the net effective downstream resistance, and the predictor corrector method predicts the net effective resistance at each centerline point (using reverse breadth first search from leaves to ostium), corrects the ostial flow (depth first search from ostium to leaves), and continues until they are consistent. For the centerline point corresponding to the outlets, the resistances are set equal to the boundary conditions, and for the rest of the centerline points, the resistances correspond to the geometric resistance [30] [31]. The flow-split correction factor, is the deviation of the assumption (flow is inversely proportional to the net effective downstream resisttance) in the original patient-specific geometry, defined at each centerline point and extracted from the patient-specific hyperemic simulations. The threshold for switching from a CFD-derived approach to a reduced order model, is chosen to be .

BUILDING PHYSICS DRIVEN RESPONSE SURFACE

Solve the HFHC model in four configurations.
for  do
     for  do
         
         
         if  then
              
              
         end if
     end for
end for

SOLVING PSROM

Initialize n = 0, ostial flow (), the superscript denoting the iteration number,
while  or  do
     n = n + 1
     
     enqueue(leaf)
     while ( do queue not empty)
         rbfs(pop queue)
     end while
     ;
end while
function dfs(current_node = m)
     if is_outlet(current_node) then return
     end if
     if is_branch(current_node) then ;
         dfs()
         dfs()
     else
         for  do .
              dfs(k)
         end for
     end if
end function
function rbfs(current_node = m)
     if is_ostium(current_node) then return
     end if
     if is_branch(current_node) then
         for  do .
              
         end for
     else
         for  do
              .
         end for
     end if
     enqueue(m-1)
end function
Algorithm 1 Overview of the algorithm to build and predict hemodynamics using a patient-specific reduced order model.

3 Results

In this section, we demonstrate the accuracy and speed of the proposed method against the HFHC computational fluid dynamics model. The detailed protocol for generating the ground truth is provided in the following section.

3.1 Validation protocol

First, we chose a cohort of patients. We identified lesions in the coronary tree of each patient, and characterized them as single (focal), ostial, bifurcation or serial lesion. Focal lesions are single isolated lesions present in a vessel path from ostium to the leaves. Bifurcation lesions go across a branch, while ostial lesions occur at the junction of aorta and coronary arteries. Three or more lesions in a vessel path were characterized as serial lesions. Any lesion type that is not characterized as one of the four listed above were not considered as candidates for the validation study.

For each patient and lesion type, we performed a virtual remodeling of the lesion and replaced the geometry within the lesion with the idealized geometry. The rest of the patient’s coronary artery tree was then blended with the remodeled region to yield a new coronary segmentation for the patient. We then performed a CFD simulation on this geometry which serves as the ground truth. The PSROM was run on the modified geometry and then compared against the ground truth simulation values.

In order to obtain sufficient statistical power, we chose thirteen hundred and forty one (1341) patients. Ten of these failed to yield a ground truth due to meshing or convergence issues. Therefore, a total of thirteen hundred and thirty one (1331) cases were available and were used to evaluate the PSROM.

The cases were gathered from centers in the US, Japan, EU, and UK and were not used for development or validation of the algorithm. Lesions that were modified had to exhibit a lumen narrowing of at least 30% using the metric defined in Sankaran et al. [29] with being used for the healthy or ideal radius. Samples of the different lesion types along with their idealized lumen geometry are shown in Figure 5.

For each case, we evaluated results at salient points downstream of the lumen modification but avoided bifurcation or pressure recovery regions. In some cases, more than one comparison location was identified. In total, 2578 comparisons across 1331 patient-specific models were performed.

In order to demonstrate equivalence, we tested that the bias is within

and the standard deviation is within

. The numbers are motivated primarily by a mesh independence study on the ground truth values, which demonstrated that if we were to replace

with a very fine mesh, the 95% confidence interval of error is within 0.03, which translates to a standard deviation of

. To demonstrate equivalence in bias, we used a two one-sided t-test (TOST) and to demonstrate equivalence in standard deviation, we used a chi-squared test. The null hypothesis for the TOST test is the absolute difference between the proposed method and ground truth is less than

. The null hypothesis for the chi-squared test is that the standard deviation is greater than 0.02.

Figure 5: Patient lumen segmentation (solid) in a section of the vessel superimposed with the ideal lumen segmentation (meshed) for (A) ostial lesion, (B) bifurcation lesion, (C) focal (single) lesion, and (D) serial lesion.

3.2 Performance

The overall bias in the evaluated dataset was 0.00524, with a corresponding p-value using TOST (two one-sided t-test) of 1e-3. The standard deviation in the evaluated dataset was 0.0147, with a corresponding p-value using chi-squared test of 1e-3. The correlation coefficient between and PSROM was 0.982 (95% CI 0.981-0.984, p0.001). The 95% Bland-Altman limits of agreement between and was (-0.020, 0.031) as seen in Figure 7. These are within the 95% confidence intervals of measurement reproducibility. Slope and intercept were observed to be 0.972 and 0.030 respectively.

The mean difference of before and after modification was 0.111 with a standard deviation of 0.121. The maximum change in before and after lumen modification was 0.70. For those locations where the change in before and after modification was greater than 0.1 (N = 1082), the bias and standard deviation were 0.006 and 0.020 respectively. For locations where the change in FFR using PSROM before and after modification was less than 0.1 (N = 1496), a bias of 0.004 and standard deviation of error of 0.009 was observed. This also demonstrates that the proposed method is able to predict for large and small changes before and after modification. Ninety five percent of all the comparisons had an error of less than 0.03 compared to .

Figure 6: (left) Run times of the proposed method with a median of 0.53 seconds, (right) A comparison of the average and standard deviation (stdev) of the full solver and the PSROM method plotted in a logarithmic scale.

The average time taken to execute the PSROM algorithm for the prediction step was sec (median 0.53 sec). The maximum time taken was 2.22 sec with a 75th percentile value of 0.75 sec. Figure 6 shows the histogram of duration of the PSROM algorithm on an iPad Pro 10.5 device. Average and standard deviation of run times across the entire population are shown for the ground truth and the PSROM methods (in logarithmic scale).

lesion type sample size bias standard deviation 95% CI
bifurcation 782 0.0051 0.0171 (-0.028, 0.039)
ostial 78 0.0089 0.0153 (-0.021, 0.039)
focal 465 0.0030 0.0149 (-0.014, 0.019)
serial 1253 0.0059 0.0084 (-0.023, 0.035)
overall 2578 0.0052 0.0148 (-0.020, 0.031)
Table 1: Performance of the proposed method categorized by lesion type. Bias, standard deviation and 95% confidence intervals are reported along with the sample size.
range N bias standard deviation
74 0.0160 0.0261
84 0.0035 0.0394
151 0.0081 0.0214
309 0.0069 0.0165
558 0.0068 0.0149
1402 0.0035 0.0080
Table 2: Performance of the proposed method stratified by . Bias and standard deviation are reported along with the sample size.

Further, the results were analyzed based on stratification by lesion type. For bifurcation lesions (N = 782), bias was 0.005 (p1e-6) and standard deviation of error was 0.017 (p1e-6). For serial lesions (N = 1253), bias was 0.006 (p1e-6) and standard deviation of error was 0.008 (p1e-6). For focal lesions (N = 465), bias was 0.003 (p1e-6) and standard deviation of error was 0.015 (p1e-6). For ostial lesions (N = 78), bias was 0.009(p=0.25) and standard deviation of error was 0.015 (p=0.004). The N corresponds to the number of comparisons, which may include more than one per patient. Table 1 summarizes these results.

Figure 7: (A) Scatterplot comparing the ground truth () solution obtained by solving the Navier-Stokes equations (HFHC model) to the solution obtained using the patient-specific reduced order model (). (B) Bland-Altman plot showing the error as a function of the average values, along with the bias and 95% limits of agreement, (C-F) Scatterplot comparing to in modifications with ostial, bifurcation, serial and single (focal) disease respectively.

Table 2 summarizes performance of the algorithm when data was stratified based on in different ranges (or buckets). Each bucket has a different sample size which is expected based on a typical randomized patient population. The bias is the highest in the bucket with smallest and the standard deviation is highest in the bucket.

Figure 8

demonstrates error properties of the FFR results of the PSROM method. It demonstrates that the error distribution is centered close to a mean of zero, and has a normal distribution around the mean. In addition, the figure also demonstrates the added value of the proposed method against a ROM either not accounting for the vessels with low flow difference or not accounting for the pressure recovery.

Figure 9 shows pressure tracing in three vessels with serial lesions with superimposed traces of the of the patient model before and after modification and . These examples demonstrate the ability of the PSROM method to capture the pressure drop in tandem lesions with higher flowrate post lumen modification.

Figure 8: Illustration of the properties of the error between and distributions. (A) Error distribution in different ranges. (B) error histogram in the evaluation dataset, (C-F) error distributions corresponding to (red) unconstrained, (orange) hybrid based on a ROM cutoff, and (green) hybrid with a pressure recovery fix term, which was used for all the results in this paper.

3.3 Validation against invasive data

We also clinically validated the PSROM method using invasive data obtained from percutaneous coronary intervention (PCI) or stenting procedures. PCI is the process of inserting a catheter with a balloon attached to restore and attain ideal lumen geometry on a section of the blood vessel. The objective of this procedure is to relieve symptoms of ischemia and mitigate risk of coronary artery disease. Some of the biggest challenges of PCI is to know a-priori which lesions have the largest impact on FFR (for serial lesions), how stenting may restore blood flow (for all lesion types, with diffuse atherosclerosis potentially not increasing coronary blood flow even after restoring vessel caliber), or how to size the stent (radius and length). The PSROM method was first applied to a patient with serial lesions where the invasive FFR was also measured [21]. Two lesions (A - proximal and B - distal) were identified in the left arterior descending (LAD) artery. Invasive revascularization of B and a combination of A and B were performed. We applied the PSROM, blinded, to these two configurations, and demonstrated excellent match for the two scenarios (lesion B - invasive 0.77, PSROM 0.79, and lesion A + B - invasive 0.85, PSROM 0.88). We also applied the PSROM method to a cohort of thirteen patients with serial lesions who underwent PCI [6]. The use of PSROM improved the correlation coefficient of the comparison to 0.75 over a value of 0.44 using conventional FFRct (p 0.001). Further, the bias and coefficient of variation of PSROM was 0.01 and 7% respectively, compared to 0.05 and 37% using pre-PCI values.

Figure 9: FFR tracing from the ostium to a distal end of a modified vessel, comparing the original trace, PSROM prediction and ground truth.

4 Discussion

We formulated and developed a method, called PSROM, that enables fast assessment of hemodynamics in response to modifying a patient-specific geometry. The method was built on the two principles. First, resistance of lumen geometry to flow was modeled to have a linear relationship with the flowrate, with the slope and intercept being functions of arterial geometry. This was primarily motivated by the steady-state solution to a 1-dimensional Navier-Stokes equation that has a linear and quadratic term for pressure loss. Second, flow split at bifurcations was modeled to be inversely related to the net downstream resistance of the respective vessels. Net downstream resistance accounts for both the geometric resistance as well as the boundary conditions. Additionally, we assumed that the idealized model represents the largest possible dilation of the CT derived lumen geometry, though the method outlined in this work is applicable even if the strategy to define idealized model is different.

We initially built a local space around a patient-specific geometry that was constrained by a maximum possible geometric expansion using a patient-specific idealized model. The idealized model was derived from the patient’s lumen segmentation by fitting a radius profile that is non-increasing from the root to the leaves. The idealized model represents a patient-specific lumen segmentation without regions of lumen narrowing that is the closest in L-1/2 norm with the original lumen segmentation. Further, we used a sampling of different microvascular properties (dilation) to derive the properties of the affine model. In this work, we define this space using four configurations, two different lumen segmentations (patient and idealized model) and two boundary conditions each. The two geometric configurations were chosen to be the extrema of the lumen segmentation (patient and idealized). The boundary conditions were chosen to be hyperemia, and another lower resistance that we call superemia, which has 40% less resistance and allows more flow through the vessels. This, in turn, enables the PSROM to parameterize the intercept and slope. A physics driven response surface was used to represent the resistance of geometry to flow within this patient-specific exploration space.

We deviated from the formulation to consider two special situations - (i) if the vessel is severely diseased such that flowrate is limited by geometry, then there is insufficient flow separation between the hyperemic and superemic conditions. In such situations, we used a model that is not derived from the HFHC simulations, but rather uses the analytical reduced order model (1D Navier-Sokes equations) constrained by the hyperemic patient and ideal model hemodynamics, (ii) in the sub-regions within the coronary tree where pressure recovers, i.e. pressure increases as we traverse the vessel tree from the proximal to distal section. These are regions where a reduced order model interpolation does not work because geometric resistance is negative. Therefore, we replace the interpolatory behavior by using an analytical model in the “pressure recovery” zone (that we define in this work as within 20 mm from the stent). The results showed that these developments were needed to achieve the desired performance.

In order to accurately capture pressure recovery, models for smaller scale fluid structures such as eddies and turbulence need to be captured. The value of 20 mm used in this paper was based on empirical observation in the development dataset, but this can be improved further by understanding the link between lumen geometry, flowrate and the resultant pressure recovery region length and extent.

The method we proposed was validated against a large cohort of patients with different lesion characteristics: namely, bifurcation, serial, single and ostial lesions. The method performed very well against the ground truth CFD simulations in over 2500 evaluation points. This was enabled by deriving the ROM directly from set of CFD simulations, thereby negating the scenarios where other low dimensional models perform poorly due to a combination of various factors such as (i) assuming a parabolic velocity profile, (ii) approximating geomety by radius, and/or (iii) failing to account for turbulent losses. Leveraging the patient-specific model to model pressure recovery, except the region just downstream of the modified region, also helped improve the performance over traditional reduced order models. The method also performed well against invasive data on a cohort of patients who underwent PCI. To our knowledge, no published method in literature has demonstrated such good performance. The technique was recently cleared by the U.S. Food and Drug Administration (FDA) for clinical use.

The proposed method showed significant discriminatory power in predicting the post modification compared to the ground truth three dimensional CFD simulations. Additionally, we demonstrated a good correlation coefficient, Bland-Altman limits of agreement, and agreement across lesion types. Two sets of clinical studies [6, 21], where the invasive measurements were blinded, showed promising results when compared with PSROM.

The method performed well in identifying “silent lesions”, those that are masked by another severe lesion in the same vessel path. Relieving the lumen narrowing in a flow-limiting lesion results in more flow through the tandem (less severe) lesion (refer Fig. 9). Consequently, this will result in a higher pressure drop across the unmodified lesion. The results we obtained demonstrated that the proposed method is able to capture this increase in pressure drop accurately.

The proposed method does not fully account for changes to the patient post-intervention such as potential plaque shift, side branch jailing, and the material used for lumen dilation. Further, the changes to the microvascular resistance are modeled only with regards to passive remodeling of the vessel wall (relationship of area to flow). Finally, the performance in modifications with a large change was not on par with the performance in healthier sections.

We believe that the method outlined in this work has the potential to aid interventional procedural planning by its ability to virtually model different scenarios before performing the actual intervention. This could potentially increase cath lab efficiency by reducing unnecessary procedures, procedural duration, total radiation dose as well as highlight scenarios that might be missed by visual inspection of an angiogram.

References

  • [1] E. Boileau and P. Nithiarasu (2015) One-dimensional modelling of the coronary circulation. application to noninvasive quantification of fractional flow reserve (FFR). Computational and Experimental Biomedical Sciences: Methods and Applications (), pp. 137–155. Cited by: §1.
  • [2] J. C. Butcher (2016) Numerical methods for ordinary differential equations. John Wiley and Sons (), pp. . Cited by: §2.6.
  • [3] P. S. Douglas, B. De Bruyne, G. Pontone, M. R. Patel, B.L. Norgaard, and R. A. B. et al. (2016) Outcomes of FFRct-guided care in patients with suspected coronary disease: the PLATFORM study. Journal of the American College of Cardiology 68 (), pp. 435–445. Cited by: §1.
  • [4] M. S. O. et al (1999) Structured tree outflow condition for blood flow in larger systemic arteries. American journal of physiology-Heart and circulatory physiology 276 (1), pp. H257–268. Cited by: §1.
  • [5] A. R. et al. (2019) Automated image acquisition for assisting a user to operate an ultrasound device. U.S. Patent Application (15/626), pp. 423. Cited by: §1.
  • [6] B. et al. (2019) Predicting the physiological effect of revascularization in serially diseased coronary arteries: clinical validation of a novel ct coronary angiography–based technique. Circulation: Cardiovascular Interventions 12 (2), pp. e007577. Cited by: §3.3, §4.
  • [7] E. B. et al. (2015) A benchmark study of numerical schemes for one‐dimensional arterial blood flow modelling. International journal for numerical methods in biomedical engineering 31 (10), pp. e02732. Cited by: §1.
  • [8] E. B. et al. (2018) Estimating the accuracy of a reduced‐order model for the calculation of fractional flow reserve (FFR). International journal for numerical methods in biomedical engineering 34 (1), pp. e2908. Cited by: §1.
  • [9] J. T. C. S. et al. (2015) Fast and accurate pressure-drop prediction in straightened atherosclerotic coronary arteries. Annals of biomedical engineering 43 (1), pp. 59–67. Cited by: §1.
  • [10] J.T.C. S. et al. (2014) Geometry-based pressure drop prediction in mildly diseased human coronary arteries. Journal of biomechanics 47 (8), pp. 1810–1815. Cited by: §1.
  • [11] K. M. et al. (2010) Atlas-based reduced models of blood flows for fast patient-specific simulations. International Workshop on Statistical Atlases and Computational Models of the Heart (Springer, Berlin, Heidelberg), pp. . Cited by: §1.
  • [12] L. I. et al. (2016) A machine-learning approach for computation of fractional flow reserve from coronary computed tomography. Journal of Applied Physiology 121 (1), pp. 42–52. Cited by: §1.
  • [13] M. .R. P. et al. (2019) 1-year impact on medical practice and clinical outcomes of FFRct: the ADVANCE registry. JACC: Cardiovascular Imaging p.2978 (), pp. . Cited by: §1.
  • [14] R. S. Driessen. et al. (2019) Head-to-head comparison of FFRct against coronary CT angiography and myocardial perfusion imaging for the diagnosis of ischemia. Journal of the American College of Cardiology 73 (2), pp. 161–173. Cited by: §1, §2.1.
  • [15] Y. O. et al. (2016) Investigation of blood flow in the external carotid artery and its branches with a new 0D peripheral model. Biomedical engineering online 15 (1), pp. 16. Cited by: §1.
  • [16] I.E. V. et. al. (2006) Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Computer methods in applied mechanics and engineering 29 (32), pp. 3776–3796. Cited by: §2.
  • [17] S. S. et. al. (2012) Patient-specific multiscale modeling of blood flow for coronary artery bypass graft surgery. Annals of biomedical engineering 40 (10), pp. 2228–2242. Cited by: §2.
  • [18] L. Formaggia, J. F. Gerbeau, F. Nobile, and A. Quarteroni (2001) On the coupling of 3D and 1D Navier–Stokes equations for flow problems in compliant vessels. Computer Methods in Applied Mechanics and Engineering 191 (6-7), pp. 561–582. Cited by: §1.
  • [19] L. Formaggia, D. Lamponi, and A. Quarteroni (2003) One-dimensional models for blood flow in arteries. Journal of engineering mathematics 47 (3-4), pp. 251–276. Cited by: §1.
  • [20] Y. Huo, M. Svendsen, J. S. Choy, Z. D. Zhang, and G. S. Kassab (2011) A validated predictive model of coronary fractional flow reserve. Journal of the Royal Society Interface 9 (71), pp. 1325–1338. Cited by: §1.
  • [21] A.R. Ihdayhid, A. White, and B. Ko (2017) Assessment of serial coronary stenoses with noninvasive computed tomography-derived fractional flow reserve and treatment planning using a novel virtual stenting application. JACC: Cardiovascular Interventions 10 (24), pp. e223–e225. Cited by: §3.3, §4.
  • [22] A. I. Khuri and S. Mukhopadhyay (2010) Response surface methodology. Wiley Interdisciplinary Reviews: Computational Statistics 2.2 (), pp. 128–149. Cited by: §1.
  • [23] H.J. Kim, I. E. Vignon-Clementel, J.S. Coogan, C.A. Figueroa, K.E. Jansen, and C.A. Taylor (2010) Patient-specific modeling of blood flow and pressure in human coronary arteries. Annals of Biomedical Engineering 38 (10), pp. 3195–3209. Cited by: §2.
  • [24] K.H. Kim, J.H. Doh, B.K. Koo, J.K. Min, A. Erglis, and H.M. Y. et al. (2014) A novel noninvasive technology for treatment planning using virtual coronary stenting and computed tomography-derived computed fractional flow reserve. JACC Cardiovascular Interventions 7 (1), pp. 72–78. Cited by: §1.
  • [25] C. Leibowitz (2017)

    Arterys: deep learning for medical imaging

    .
    Demystifying Big Data and Machine Learning for Healthcare (), pp. 169–173. Cited by: §1.
  • [26] M. D. McKay, R. J. Beckman, and W.J. Conover (1979) Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), pp. 239–245. Cited by: §2.3.
  • [27] V. Milišić and A. Quarteroni (2004) Analysis of lumped parameter models for blood flow simulations and their relation with 1D models. ESAIM: Mathematical modelling and numerical analysis 38 (4), pp. 613–632. Cited by: §1.
  • [28] B.L. Norgaard, J. Leipsic, S. Gaur, S. Seneviratne, B. S. Ko, H. Ito, J.M. Jensen, L. Mauri, B. De Bruyne, H. Bezerra, K. Osawa, M. Marwan, C. Naber, A. Erglis, S. J. Park, E.H. Christiansen, A. Kaltoft, J. F. Lassen, H. E. Btker, and S. Achenbach (2014) Diagnostic performance of non-invasive fractional flow reserve derived from coronary CT angiography in suspected coronary artery disease: the NXT trial. Journal of American College of Cardiology 63 (12), pp. 1145–1155. Cited by: §1, §2.1.
  • [29] S. Sankaran, L. Grady, and C.A. Taylor (2015) Fast computation of hemodynamic sensitivity to lumen segmentation uncertainty. IEEE Transactions on Medical Imaging 34 (12), pp. 2562–2571. Cited by: §3.1.
  • [30] S. Sankaran, L. Grady, and C.A. Taylor (2015) Impact of geometric uncertainty on hemodynamic simulations using machine learning. Computer Methods in Applied Mechanics and Engineering 297 (), pp. 167–190. Cited by: §2.5, §2.5, §2.8.
  • [31] S. Sankaran, H.J. Kim, G. Choi, and C.A. Taylor (2016) Uncertainty quantification in coronary blood flow simulations: impact of geometry, boundary conditions and blood viscosity. Journal of biomechanics 49 (12), pp. 2540–2547. Cited by: §2.5, §2.8.
  • [32] S. Sankaran and A.L. Marsden (2010) The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow. Physics of Fluids 22 (), pp. 121902. Cited by: §1.
  • [33] S. Sankaran and A.L. Marsden (2011) A stochastic collocation method for uncertainty quantification and propagation in cardiovascular simulations. Journal of Biomechanical Engineering 133, pp. 031001–1. Cited by: §2.3.
  • [34] B. N. Steele, D. Valdez-Jasso, M. A. Haider, and M. S. Olufsen (2011) Predicting arterial flow and pressure dynamics using a 1D fluid dynamics model with a viscoelastic wall. SIAM Journal on Applied Mathematics 71 (4), pp. 1123–1143. Cited by: §1.
  • [35] C. A. Taylor, T. J. R. Hughes, and Zarins,C. K. (1998) Finite element modeling of blood flow in arteries. Computer Methods in Applied Mechanics and Engineering 158 (1-2), pp. 155–196. Cited by: §2.
  • [36] C.A. Taylor, M.T. Draney, J.P. Ku, D. Parker, B.N. Steele, K. Wang, and C.K. Zarins (1999) Predictive medicine: computational techniques in therapeutic decision-making. Computer aided surgery 4 (5), pp. 231–247. Cited by: §1.
  • [37] C.A. Taylor, T. A. Fonte, and J. K. Min (2013) Computational fluid dynamics applied to cardiac computed tomography for noninvasive quantification of fractional flow reserve: scientific basis. Journal of American College of Cardiology 61 (22), pp. 2233–2241. Cited by: §1, §1, §2.
  • [38] D. Xiu and J. S. Hesthaven (2005) High order collocation methods for the differential equation with random inputs. Journal of Scientific Computing 27 (), pp. 1118–1139. Cited by: §2.3.