Patient-Specific Characterization of Breast Cancer Hemodynamics Using Image-Guided Computational Fluid Dynamics

1). Segmentation and vessel detection

The registered, high spatial resolution DCE-MRI data were used for morphological analysis using a framework previously developed. While details are presented elsewhere [39], we briefly summarize the salient points. The lesions were segmented from the surrounding breast tissue using a fuzzy c-means clustering method [40]. (The same method was also adapted for segmentation of fibroglandular and adipose tissues.) Before vessel segmentation, histogram normalization and enhancement were performed by applying a local-statistics-based, intensity transfer function to the subtraction of the pre-from the post-contrast, high-spatial resolution images [41]. The 3D vasculature was then segmented from the surrounding tissue by applying a Hessian filter to the enhanced images [42], from which a probability for each voxel to belong to the vasculature was constructed. Finally, a lowest-cost tracking algorithm automatically detected vessels that have a high probability of contacting the lesion, as well as the path yielding the most likely physical connection between the vessel and tumor [39].

2). Pharmacokinetic modeling

The ultrafast DCE-MRI data were analyzed by a Patlak model [43] modified with a delayed bolus-arrival time (BAT):

where C_t (t) and C_p(t) are the time courses of the concentration of contrast agent in the tissue and the local plasma, respectively, and C_p,pop(t) is the population arterial input function measured from the internal thoracic arteries. Thus, the BAT, K^trans (volume transfer coefficient, representing local vascular perfusion and permeability), and v_p (plasma volume fraction) can be estimated for each voxel. More detailed explanation of the identification of C_p(t) and pharmacokinetic modeling can be found in [39]. Examples of DCE time courses and fitting can be found in the supplemental materials.

DW-MRI data were analyzed by standard methods [44] to yield parametric maps of the apparent diffusion coefficient (ADC). As the DW-MRI data were acquired with suppression of the fat signal (as is commonly done for DW-MRI of the breast), directly mapping the ADC values of that tissue component is challenging. Thus, for adipose tissue, a literature value (ADC = 1.91×10⁻⁵ mm²/s) was assigned to all voxels labeled as adipose tissue after segmentation [45]. In addition, a small number of unphysical ADC values (i.e., ADC values below 0.0 or above 3.0×10⁻³ mm²/s) may also exist in fibroglandular tissue due to signal noise. Thus, any voxels within the fibroglandular or tumor tissues displaying unphysical values were replaced by the average of their valid nearest-neighbor (i.e., 26 voxels in 3D). Finally, Gaussian smoothing with a standard deviation of 1.0 was applied at the edges of tissue regions to reduce unphysical discontinuities.

3). Vessel scaling and mesh generation

Numerical meshes of the whole breast including the lesion and vasculature were generated from segmented masks obtained from section II.C.1). Due to the relatively small radius of breast vessels compared to the imaging resolution, the directly segmented vessel mask suffers from significant partial volume effects and cannot accurately represent the size of vessels. Thus, an adaptive vessel scaling algorithm was designed by integrating information from the segmented vessel mask and the v_p map.

First, the segmented vascular mask is skeletonized to its centerline (Fig. 2(a)). We identify the branching points, terminal ends, and vessel segments of the vascular network, and calculate orientations of each vessel. A gap filling process is applied with a similar strategy as the tumor-associated vessel tracking, so that the vascular network becomes a single connected graph. Furthermore, a moving average of length seven is applied to smooth the vessel centerlines. Second, the mask and v_p map are up-sampled to an isotropic resolution of 200 microns (Fig. 2(b)). Third, at each specific location along the vessel centerline, a neighborhood in the up-sampled grid is identified using the vascular orientation and a step width dl of half of the up-sampled voxel size. This neighborhood represents the region in the original vessel mask covering the section of vessel at this location (Fig. 2(c)). We consider this section to be a cylinder with height dl and an unknown radius R, while its volume should be approximately equal to the v_p value for this neighborhood multiplied by voxel volume in the refined grid:

where v_p(i) is the measured plasma volume fraction for the i^th voxel in this identified neighborhood, and V_vox(i) is the volume of the corresponding voxel, which is a constant 0.008 mm³.

Eq. (18) enables the calculation of the local vascular radius. The voxels within the range of the cylinder with height dl and calculated radius R are considered the vascular region (Fig. 2(d)). The algorithm just described is applied stepwise to the entire vascular network to properly scale the vessel mask throughout the domain (Figs. 2(e) and 2(f)). Finally, a Matlab toolkit, iso2mesh [46], is used to convert the binary mask to surface and volumetric numerical meshes of the vasculature. The same toolkit is also used for generating meshes of the lesion and whole breast regions from their respective masks. The extravascular tissue mesh is then generated by applying the Boolean difference operation (‘surfboolean’ function in iso2mesh) between the breast contour and the vascular surfaces.

4). Parameter scaling

Parametric maps of ADC and K^trans inform spatial-resolved (and patient-specific) values of tissue hydraulic conductivity, κ, and vascular permeability, L_p, respectively. We assume a simple linear scaling where the median of the scaled MRI parameter matches the median of the corresponding material property value from the literature.

5). Boundary conditions and direction of blood flow

Appropriate boundary conditions are necessary on both the boundaries of 3D tissue domain, ∂Ω_t, as well as the terminal ends of pseudo-1D vascular network. Three subdomains of the tissue boundary are identified: 1) the surface contour of the breast, ∂Ω_t,1, 2) the chest wall, ∂Ω_t,2, and 3) the exterior surface of the vascular domain, ∂Ω_t,3 = ∂Ω_t∩∂Ω_v. In this study, Dirichlet boundary condition for pressure (i.e., constant pressure p_t = 2×10⁴ g·cm⁻¹·s⁻²) [31] is assigned on ∂Ω_t,1 and ∂Ω_t,2, and Neumann condition is assigned on ∂Ω_t,3 (i.e., Eq. (14)).

The Dirichlet conditions on the pressure field at the terminal ends of the vascular network, as well as the direction of blood flow, are automatically determined using an optimization procedure modified from the work of Fry et al. [47]. Briefly, blood flow is first modeled with a simplifying assumption (i.e., non-permeable vessels and a uniform vascular radius) and randomly assigned flow directions. Because the boundary condition on the terminal ends is initially unknown, the simulation is implemented as an optimization process with flux continuity constraints at the vascular junction points. The optimization goal is to approach the target pressures at each node and the target shear stresses in each vessel segment. Thus, we minimize the object function L in Eq. (19):

with

where k_p and k_τ weight the relative contribution of the pressure and shear stress terms, respectively; p_v,k is the pressure at the k^th of N total vascular nodes; δ_ik is the Kronecker delta; l_j and τ_j are the length of and shear stress in the j^th vessel segment, respectively, of S total vessels; S_k the set of vessels connected to k^th node; and S_ik the set of vessels connecting the i^th and k^th nodes. The third term in Eq. (19) arises through the constraint of flux continuity at the set of internal nodes, which are enforced with a standard Lagrange multiplier approach, where λ_i represents the Lagrange multiplier for the i^th node. The target value of wall shear stress is set to be ∣τ₀∣ = 15 g·cm⁻¹·s⁻², and target pressure p_v,0 = 4.13 × 10⁴ g·cm⁻¹·s⁻² [48]. The viscosity of blood, μ, is set to be 0.04 g·cm⁻¹·s⁻¹ [31][49].

After each iteration of this simulation, the bolus arrival-time, BAT^(estimate), is then estimated and compared with the measured bolus arrival-time map, BAT^(measure). Finally, simulated annealing (starting from a small k_τ and doubling the ratio k_τ/k_p during each iteration to avoid biases toward the initial randomly assigned directions) is applied to adjust the assignment of flow directions until the normalized difference between BAT^(estimate) and BAT^(measure), M, is minimized; i.e., we minimize:

This optimization procedure determines the flow direction and pressure at all terminal ends.

6). Numerical implementation

To numerically solve the 1D-3D coupled simulation system, (i.e., Eqs. (10) - (14)), an iterative algorithm which combines the finite difference method (FDM) and the finite element method (FEM) was implemented using the Python FEniCS library [50] as follows (definitions can be found in section II.A):

Step 1: Initialize pressure at exterior surface of vessel, p_t,ev ⁽⁰⁾.

Step 2: Solve Eq. (10) for p_v⁽¹⁾ with the FDM on the vessels,

where for each j ∈ {junction points},

Step 3: With the estimated p_v⁽¹⁾ and Eq. (14), solve Eq. (13) for p_t⁽¹⁾ using the FEM with tetrahedral elements and a test space of second-order, continuous, piecewise polynomials:

Step 4: Update p_t,ev⁽¹⁾ by averaging the estimated p_t⁽¹⁾ along the perimeter of vessel at position l.

Step 5: Repeat Steps 2 – 4, until p_v⁽ⁿ⁾, p_t⁽ⁿ⁾ converge according

Once this system (i.e., Eqs. (10) - (14)) is solved, it provides the fields of steady-state blood pressure, p_v, and interstitial pressure, p_t, throughout the domain. The local blood flow rate, Q_v, and interstitial flow velocity, u_t, can now be calculated based on Eqs. (1) and (2), respectively.

7). Computation time

The CFD simulation of each patient case takes approximately three minutes to converge on our lab’s server (CPU: 4× Intel Xeon E5-4627 v4 2.6GHz, 25M Cache, 8.0 GT/s; Memory: 256 GB), and approximately ten minutes on a standard laptop (MacBook Pro with CPU: 2.7GHz Intel Core i5, Memory: 8GB 1867MHz DDR3), even without using parallel solvers. The pre-processing are actually the major time-consuming steps, they include: calculating the T₁ map, pharmacokinetic fitting, meshing of the 3D tissue domain, and defining the coupling relationship between 1D and 3D location coordinates. This takes approximately two hours for each patient with our current implementation. It should be noted, though, that these pre-processing steps are voxel- or node-wise computations, which can be directly parallelized. Further discussion on the robustness and uncertainty of the modeling pipeline, as well an assessment of error propagation from the input parameters to the computation outputs can be found in the supplemental materials.

1). Image processing results

Fig. 3 shows the direct outputs of image processing for one example patient with a malignant lesion (panels (a) – (e)) and one example patient with a benign lesion (panels (f) – (j)). On the left side, panels (a) and (f) show 3D reconstructions of the vascular skeleton for the whole breast, lesions, and detected tumor-associated vessels. For the malignant lesion, seven tumor-associated vessels were detected, with one inlet and six outlets. The benign lesion has three tumor-associated vessels detected, with one inlet and two outlets. On the right side of the figure are the measured BAT (panels (b) and (g)), v_p (panels (c) and (h)), K^trans (panels (d) and (i)) and ADC (panels (e) and (j)) parametric maps in diseased breasts. The median of parameters for the example malignant tumor are BAT = 9.98 s, v_p = 1.45×10⁻², K^trans = 0.34 min⁻¹, and ADC = 1.15×10⁻³ mm² s⁻¹; parameters for the benign tumor are BAT = 11.10 s, v_p = 1.00×10⁻⁴, K^trans = 0.04 min⁻¹ and ADC = 9.60×10⁻⁴ mm² s⁻¹.

2). Mesh generation and material properties assignment

Fig. 4 presents the results of mesh generation and properties assignment for the example malignant (panels (a) – (d)) and the benign (panels (e) – (h)) cases respectively. Panels (a) and (e) visualize the numerical meshes of vasculature and the tumor regions. K^trans maps are scaled to L_p and assigned to the vessel mesh (panels (c) and (g)). The literature values of L_p fall between 2.5×10⁻¹⁰ ~ 2.1×10⁻⁹ g⁻¹ cm² s with a median of 1.2×10⁻⁹ g⁻¹ cm² s [29][51]. The median of K^trans within the breast vasculature for all patients is 0.16 min⁻¹. Thus, a scalar coefficient of k₁ = 7.6 × 10⁻⁹ is applied to scale the K^trans values. Similarly, ADC maps are scaled to κ and assigned to the tissue mesh (panels (d) and (h)). The literature values of κ are between 10⁻¹² ~ 10⁻¹⁰ g⁻¹ cm³ s yielding a median value of 1.0 × 10⁻¹¹ g⁻¹ cm³ s [29][52]. The median of ADC for the whole breast tissue for all patients is 1.91×10⁻⁵ mm² s⁻¹. A scalar coefficient of k₂ = 5.0 × 10⁻⁷ is applied to the ADC values.

3). Automatic determination of the direction of blood flow

Fig. 5 illustrates the results of the hierarchical optimization procedure for determining the boundary condition at vascular terminal ends and the blood flow direction in the example patient with benign lesion. The objective function M reaches a minimal value of 0.61 at the 63^th iteration of the optimization procedure (panel (a)). The resulting distribution of estimated BAT shows a DICE similarity of 0.83 with the measured BAT (panel (b)). Spatial-resolved BATs are visualized in panel (c).

1). Whole breast pressure fields

Fig. 6 presents solutions of the image-based flow simulation for the diseased breasts of the two example patients. Panels (a) and (e) show the blood pressure along the vasculature network for the malignant and benign diseased breasts, respectively, where major inlet arteries (indicated by white arrows) from the posterior and lateral sections of the breasts can be identified via the color-coding. Panels (b) and (f) show the interstitial pressure fields in the breasts. Due to the significant efflux from vessel to breast tissue, regions close to the exterior surface of the major arterial inlets manifest high interstitial pressure, with the median and range of 2.01×10⁴ (2.00×10⁴ − 2.02×10⁴) g·cm⁻¹·s⁻²; while interstitial pressure close to the venous surface is generally low, with the median and range of 1.999×10⁴ (1.995×10⁴ − 2.000×10⁴) g·cm⁻¹·s⁻², where influx from tissue back to vessels presents.

2). Tumor interstitial pressure and flow velocity fields

Enlarged views of the interstitial pressure and flow velocity fields of the tumor are also shown in Fig. 6. Panels (c) and (g) present the pressure fields in the malignant and benign tumors, respectively. Comparing these two panels, the malignant lesion tends to have higher and wider range of interstitial pressure value; the median of p_t for the malignant and benign tumors are 2.029×10⁴ g·cm⁻¹·s⁻² and 2.015×10⁴ g·cm⁻¹·s⁻², respectively. Panels (d) and (h) present the velocity vector fields in the malignant and benign tumors, with the arrow sizes scaled by the local vector magnitude. The magnitude of interstitial flow velocity in the benign tumor is smaller than that in the malignancy, where the median of ∣u_t∣ for presented malignant and benign tumors are 7.72×10⁻⁸ cm·s⁻¹ and 1.11×10⁻⁸ cm·s⁻¹ respectively. In Fig. 6(d), obvious inward interstitial flow is indicated by the velocity vectors at the superior regions of the malignant tumor, where a major inlet artery contacts the lesion. Outward flow vectors are also observed at the inferior and anterior regions of the tumor, while flow is weak in the medial and posterior regions. Fig. 6(h) shows the interstitial flow in the benign tumor. (For movies showing the interstitial pressure and flow fields, please see the supplemental materials.)

1). Spearman’s correlation between parameters

With Spearman’s correlation test performed, only the range of tumor interstitial flow velocity (range ∣u_t∣) is strongly correlated with other parameters. Specifically, ρ = 0.90, 0.92, and 0.95 for the correlations between the “range ∣u_t∣” and “median p_t”, “range q_e” and “range p_t”, respectively. Thus, the parameter “range ∣u_t∣” is excluded from further statistical tests.

2). Wilcoxon test of each hemodynamic characteristic

As summarized in Table II, the median of tumor interstitial flow velocity (median ∣u_t∣), range of blood pressure (range p_v), and extraction rate (range q_e) along the tumor-associated vessels, all have significant differences between malignant and benign lesions, based on the Wilcoxon test (P = 0.028, 0.016 and 0.040, respectively; see Fig. 7).

3). ROC analysis to evaluate diagnostic performance

We use a multivariate regression model to assess the ability of the hemodynamic characteristics to separate malignant from benign disease. Among all possible combinations of the parameters, the model generated using a combination of {median of q_e, median of ∣u_t∣ and range of p_v} performs the best. It achieves a sensitivity, specificity, and AUC of 1.0 in this proof-of-principle analysis.