1 Introduction
Patientspecific 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 patientspecific 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 timeaveraged pressure downstream of a coronary artery stenosis to a timeaveraged 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 patientspecific 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 patientspecific noninvasive 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 tubelike stenosis model, however calculation of some of the higher order geometric quantities (higher order gradients of radius or area) accurately on patientspecific 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 patientspecific 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 highfidelity 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 preset 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 preset 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 realtime 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 onedimensional NavierStokes 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 highfidelity model for the configurations where highfidelity simulations were performed.
We demonstrate that we can achieve excellent accuracy of the proposed method compared to threedimensional 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 patientspecific reduced order model. We describe the HFHC model and demonstrate how we can build a realtime 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 patientspecific anatomic geometries that are reconstructed from coronary CT scans. The NavierStokes equations are used to model blood flow and lumped parameter boundary conditions are prescribed that account for the flowpressure relationship in the microvasculature (vessels that are not visible in the CCTA and therefore not modeled). VignonClementel et al. [16] developed methods for coupling boundary conditions using a circuit analogy of the microvasculature with the NavierStokes 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 noninvasive 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 nonlinear differential operator modeling the NavierStokes 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 onedimensional 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) 
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 MonteCarlo 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 onedimensional, steady state NavierStokes 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 patientspecific 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 patientspecific geometry. The reduced order operator contains the affine resistance model that relates flowrates to pressures and the flowsplit 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.
2.4 Patientspecific ideal geometries
One of the novelties of our work is to build a patientspecific 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 nonincreasing 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 pointwise 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 patientspecific 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.
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 patientspecific 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 crosssection 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 patientspecific 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 1dimensional modeling of the NavierStokes equations [30]. This modification can be done for both the patientspecific 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 patientspecific 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 Patientspecific Reduced order model
The patientspecific 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 patientspecific 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 preset 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 patientspecific 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 nonlinear 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 patientspecific model.
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 patientspecific 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 flowsplit correction factor, is the deviation of the assumption (flow is inversely proportional to the net effective downstream resisttance) in the original patientspecific geometry, defined at each centerline point and extracted from the patientspecific hyperemic simulations. The threshold for switching from a CFDderived approach to a reduced order model, is chosen to be .
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 patientspecific 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 replacewith 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 onesided ttest (TOST) and to demonstrate equivalence in standard deviation, we used a chisquared 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 chisquared test is that the standard deviation is greater than 0.02.3.2 Performance
The overall bias in the evaluated dataset was 0.00524, with a corresponding pvalue using TOST (two onesided ttest) of 1e3. The standard deviation in the evaluated dataset was 0.0147, with a corresponding pvalue using chisquared test of 1e3. The correlation coefficient between and PSROM was 0.982 (95% CI 0.9810.984, p0.001). The 95% BlandAltman 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 .
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) 
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 
Further, the results were analyzed based on stratification by lesion type. For bifurcation lesions (N = 782), bias was 0.005 (p1e6) and standard deviation of error was 0.017 (p1e6). For serial lesions (N = 1253), bias was 0.006 (p1e6) and standard deviation of error was 0.008 (p1e6). For focal lesions (N = 465), bias was 0.003 (p1e6) and standard deviation of error was 0.015 (p1e6). 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.
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.
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 apriori 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 prePCI values.
4 Discussion
We formulated and developed a method, called PSROM, that enables fast assessment of hemodynamics in response to modifying a patientspecific 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 steadystate solution to a 1dimensional NavierStokes 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 patientspecific geometry that was constrained by a maximum possible geometric expansion using a patientspecific idealized model. The idealized model was derived from the patient’s lumen segmentation by fitting a radius profile that is nonincreasing from the root to the leaves. The idealized model represents a patientspecific lumen segmentation without regions of lumen narrowing that is the closest in L1/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 patientspecific 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 NavierSokes equations) constrained by the hyperemic patient and ideal model hemodynamics, (ii) in the subregions 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 patientspecific 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, BlandAltman 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 flowlimiting 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 postintervention 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] (2015) Onedimensional 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] (2016) Numerical methods for ordinary differential equations. John Wiley and Sons (), pp. . Cited by: §2.6.
 [3] (2016) Outcomes of FFRctguided 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] (1999) Structured tree outflow condition for blood flow in larger systemic arteries. American journal of physiologyHeart and circulatory physiology 276 (1), pp. H257–268. Cited by: §1.
 [5] (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] (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] (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] (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] (2015) Fast and accurate pressuredrop prediction in straightened atherosclerotic coronary arteries. Annals of biomedical engineering 43 (1), pp. 59–67. Cited by: §1.
 [10] (2014) Geometrybased pressure drop prediction in mildly diseased human coronary arteries. Journal of biomechanics 47 (8), pp. 1810–1815. Cited by: §1.
 [11] (2010) Atlasbased reduced models of blood flows for fast patientspecific simulations. International Workshop on Statistical Atlases and Computational Models of the Heart (Springer, Berlin, Heidelberg), pp. . Cited by: §1.
 [12] (2016) A machinelearning approach for computation of fractional flow reserve from coronary computed tomography. Journal of Applied Physiology 121 (1), pp. 42–52. Cited by: §1.
 [13] (2019) 1year impact on medical practice and clinical outcomes of FFRct: the ADVANCE registry. JACC: Cardiovascular Imaging p.2978 (), pp. . Cited by: §1.
 [14] (2019) Headtohead 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] (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] (2006) Outflow boundary conditions for threedimensional 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] (2012) Patientspecific multiscale modeling of blood flow for coronary artery bypass graft surgery. Annals of biomedical engineering 40 (10), pp. 2228–2242. Cited by: §2.
 [18] (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 (67), pp. 561–582. Cited by: §1.
 [19] (2003) Onedimensional models for blood flow in arteries. Journal of engineering mathematics 47 (34), pp. 251–276. Cited by: §1.
 [20] (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] (2017) Assessment of serial coronary stenoses with noninvasive computed tomographyderived 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] (2010) Response surface methodology. Wiley Interdisciplinary Reviews: Computational Statistics 2.2 (), pp. 128–149. Cited by: §1.
 [23] (2010) Patientspecific modeling of blood flow and pressure in human coronary arteries. Annals of Biomedical Engineering 38 (10), pp. 3195–3209. Cited by: §2.
 [24] (2014) A novel noninvasive technology for treatment planning using virtual coronary stenting and computed tomographyderived computed fractional flow reserve. JACC Cardiovascular Interventions 7 (1), pp. 72–78. Cited by: §1.

[25]
(2017)
Arterys: deep learning for medical imaging
. Demystifying Big Data and Machine Learning for Healthcare (), pp. 169–173. Cited by: §1.  [26] (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] (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] (2014) Diagnostic performance of noninvasive 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] (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] (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] (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] (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] (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] (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] (1998) Finite element modeling of blood flow in arteries. Computer Methods in Applied Mechanics and Engineering 158 (12), pp. 155–196. Cited by: §2.
 [36] (1999) Predictive medicine: computational techniques in therapeutic decisionmaking. Computer aided surgery 4 (5), pp. 231–247. Cited by: §1.
 [37] (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] (2005) High order collocation methods for the differential equation with random inputs. Journal of Scientific Computing 27 (), pp. 1118–1139. Cited by: §2.3.
Comments
There are no comments yet.