Computer simulations suggest that prostate enlargement due to benign prostatic hyperplasia mechanically impedes prostate cancer growth

Patient Data.

Anonymized patient data were obtained from the public repository for PCa imaging data at the Initiative for Collaborative Computer Vision Benchmarking website (i2cvb.github.io/) (58). Institutional review board approval and informed consent were not required for this study. We used the multiparametric MR image dataset obtained with a 3.0 Tesla Siemens Magnetom Trio TIM scanner. Patient cohort description and details on data acquisition have been previously described (59). For each biopsy-confirmed patient, this database includes T2-weighted MR images; dynamic contrast enhanced MR images; diffusion-weighted MR images; MR spectroscopic images; apparent diffusion coefficient maps; and the segmentations of the prostate, the PZ, the CG, and the tumor by an experienced radiologist. For this research, we selected a patient aged 54 y at MR date who had a large prostate of 52.81 cc harboring a localized tumor in the PZ.

Mathematical Model.

Our modeling philosophy aims at developing tumor growth models based on key phenomena and featuring a limited number of representative parameters (30). This approach holds the potential for precise heterogeneous parameterization using available longitudinal clinical and imaging data from patients and has proved successful in reproducing and predicting various types of tumor growth (30, 32⇓–34). As a preliminary step in this direction, here we developed a mechanically coupled model for patient-specific, organ-confined PCa growth and we leveraged it to perform a study on the mechanical interaction between PCa and BPH using literature data for parameter calibration.

Our model is described by the following equations:∂ϕ∂t=MDϕΔϕ−1τdF(ϕ)dϕ+χs−Aϕ[1]∂s∂t=DsΔs+S−δϕ−γss[2]∂p∂t=DpΔp+αh(1−ϕ)+αcϕ−γpp[3]∇⋅σ=0.[4]Eqs. 1–3 come from our previous model of organ-confined PCa growth (39, 40). Eq. 1 describes tumor dynamics using the phase-field method (41). The order parameter ϕ takes values from 0 in healthy tissue to 1 within tumoral regions, showing a thin diffuse interface in between. Fϕ=16ϕ21−ϕ2 is a double-well potential, which enables the stable coexistence of healthy and cancerous tissue in our model. The last two terms in Eq. 1 describe nutrient-driven growth and apoptosis (i.e., programed cell death), respectively. We consider that tumor growth is driven by a generic nutrient s that follows reaction–diffusion dynamics in Eq. 2. The reactive terms in Eq. 2 are the nutrient supply, the consumption of nutrient by the tumor, and a natural decay, respectively. Eq. 3 describes the reaction–diffusion dynamics of tissue PSA p, defined as the serum PSA concentration leaked to the bloodstream per unit volume of prostatic tissue (39). In Eq. 3, we consider that healthy tissue and cancerous tissue produce p at rates αh and αc, respectively, and the last term is a natural decay.

The mechanical stress fields generated by a growing tumor have been shown to slow down its dynamics (23⇓⇓⇓⇓–28). Previous clinically oriented, mechanically coupled models of tumor growth have incorporated this phenomenon by means of a coefficient that exponentially decreases cell mobility or net cell proliferation as mechanical stress fields intensify (30, 32⇓–34). Following a similar strategy, in Eq. 1 we modeled the mechanically induced inhibition of tumor growth throughout coefficient M, which we define in Eq. 5. Mechanical stress can be decomposed as the sum of hydrostatic stress σh (defined in Eq. 6), which tends to change the volume of the stressed body, and deviatoric stress, which tends to distort it (46). The Von Mises stress σv (46) (defined in Eq. 7) is a good measure of the distortional strain energy around the tumor that has been widely adopted in mechanically coupled models (32⇓–34). However, σv is insensitive to hydrostatic stress, which better describes the natural stress state within a growing region of tissue, such as a tumor or the CG developing BPH (23, 24, 26). Hydrostatic stress is also bound to increase globally during these growth processes because the prostate is a confined organ (1). Therefore, we propose to combine σv and the hydrostatic stress σh (46) in the definition of M following a similar approach to that of multiaxial stress-based failure criteria,M=e−β1(σv+β2σh),[5]whereσh=13σ:I=13σ11+σ22+σ33,[6]σv=σ112+σ222+σ332−σ11σ22−σ22σ33−σ33σ11+3σ122+σ232+σ1321/2,[7]where σij with i,j=1,2,3 are the components of the stress tensor σ, I is the second-order identity tensor, and β1 and β2 are constants that were adjusted to match the experimental and clinical results of previous studies on tumor growth (23⇓⇓⇓⇓–28, 32⇓–34). We use the absolute value of the hydrostatic stress because both growth-induced compressive and tensile solid stresses have been shown to impede tumor dynamics at tissue scale in vivo (23, 24).

We assume that PCa and BPH show sufficiently slow rates (1) to neglect the inertial effects in the deformation of prostatic tissue. Hence, Eq. 4 describes mechanical equilibrium as a quasistatic process. Linear elasticity has been widely accepted to model mechanical equilibrium in living tissue subjected to slow processes over short time scales (t∼1 y) (32⇓⇓⇓⇓⇓–38, 46). Hence, this paradigm provides a simple mechanical framework to compute the stresses in Eq. 5 and perform our qualitative study on the mechanical influence of BPH on PCa. The prostate is an histologically heterogeneous organ: While the CG has a higher and denser stromal component, the PZ has more abundant glandular elements with sparsely interwoven smooth muscle (1, 5, 6). Benign enlargement of the prostate tends to make the CG denser and more compact (1, 7). Therefore, the CG is normally stiffer than the PZ in patients with PCa and/or BPH (42⇓⇓–45). The growing tumor induces an internal compressive hydrostatic stress and exerts outward forces acting on the tumor border (22, 24), so we modeled the tumor mass effect asptumor=−κϕ.[8]This assumes that κ is the magnitude of a constant compressive pressure, which is admissible over short simulation times (32⇓–34, 36, 37). Volumetric growth due to BPH is usually modeled with an exponential function, but the slow growth rates justify a linear approximation over periods Δt ∼1–10 y (51⇓–53). Hence, we modeled BPH as a linearly growing pressure acting only on the CG,pBPH=−KϱgΔtVMRIHCGx,[9]where K is the bulk modulus, g is the estimated linear volumetric growth rate of the prostate due to BPH at MR date, VMRI is the volume of the prostate as measured on T2-weighted MRI, Δt is the time passed since the date of MRI, and HCGx is a Heaviside function with value 1 in the CG and 0 elsewhere. The parameter ϱ adjusts the value of the pressure pBPH to obtain the estimated CG growth rate. Taking all of the considerations above, we model prostatic tissue as a linear elastic, heterogeneous, isotropic material whose constitutive equation is given byσ=λ∇⋅uI+2μ∇su−κϕI−KϱgΔtVMRIHCGxI,[10]where σ is the stress tensor, λ and μ are the Lamé constants, u is the displacement vector, and x is the position vector. The segmentation of CG and PZ was extracted from our patient’s imaging data and mapped over the quadrature points to define heterogeneous material properties (Fig. 1).

We computed tumor volume Vϕ and serum PSA P as (39)Vϕ=∫ΩϕdΩP=∫ΩpdΩ,[11]where Ω is the prostate segmented on the patient’s MR images.

Because we are focusing on localized PCa, we imposed zero-valued Dirichlet conditions for ϕ all over the prostate boundary ∂Ω. We set natural boundary conditions for s and p. The confinement of the prostate within the pelvic area (1) was modeled with Winkler-inspired boundary conditions on the external surface of the prostate, while free displacement was enabled along the urethra; i.e.,σn=−kwuin ∂Ωexteriorσn=0in ∂Ωurethra,[12]where n is the outer normal vector to ∂Ω and kw is constant.

The initial condition of the phase field, ϕ0, was an L2 projection of the tumor segmentation extracted from the patient’s T2-weighted MR images and mapped over quadrature points or was artificially modeled with an L2-projected hyperbolic tangent field. The initial conditions for s and p were approximated with linear functions based on ϕ0 (40). Linear elasticity allows us to apply the principle of superposition and set zero-valued initial conditions for the displacements; i.e., u0=0. This means that we are computing only the displacements produced since the MR date.

To study how BPH mechanically influences PCa growth, we estimated the stress state σ0 caused by the history of BPH before the detection of PCa at MR date and leveraged it as a prestress state to compute M in Eq. 1. Hence, for this purpose, σ=σ0+σ1, where σ1 are the stresses developed since the detection of PCa at MR date. To estimate σ0, we assumed that the volume of the patient’s prostate was 20 cc at age 40 y (1, 3) and leveraged Eq. 4 with a negative value for g and Δt=14 y to approximate the undeformed, healthy state of our patient’s prostate according to the standard anatomical features (1). Because we are using linear elasticity, it suffices to reverse the sign of the obtained displacements to yield u0 and then approximate σ0 asσ0=λ∇⋅u0I+2μ∇su0−KϱgΔt0VMRIHCGxI,[13]where g is now positive and Δt0=14 y.

Parameter selection for Eqs. 1–3 has been previously discussed (39) and it is provided in Table 1. We set the values of the Young modulus of PZ and CG to EPZ=3 kPa and ECG=6 kPa, respectively (42⇓⇓–45). Because living soft tissues have a high content of water (34, 37, 42), we set the Poisson coefficients νPZ=νCG=0.40. We empirically calibrated β1=0.80 1/kPa and β2=1.50 in agreement with the experimental and clinical observations reported in previous studies on tumor growth (23⇓⇓⇓⇓–28, 32⇓–34). We selected κ=2.50 kPa to produce displacements in the order of 1 mm. The analysis of BPH and PCa deformation and the computation of the stress fields σ0 rendered g=2.34 cc/y, ϱ=2.78, and kw=0.23 kPa/mm. These values for g and ϱ reasonably agree with previous observations in the literature for our patient characteristics (1, 3, 50⇓–53).

Finally, we acknowledge that linear elasticity is an acceptable simplification for studies featuring small strains and rotations and that the strain values in the simulations shown in Figs. 3–5 are somewhat outside the admissible range for linear elasticity in a few localized regions along the tumor interface. Additionally, the prestress σ0 is just a gross estimate of the stress state in the prostate due to the history of BPH previous to PCa detection at MR. However, we are using the stresses only to compute the mechanotransductive factor to adjust tumor dynamics. This study also features other major sources of uncertainty beyond linear elasticity, such as prostate and tumor segmentation, the PCa model itself, or the mechanical boundary conditions. Still, our predictive simulations still qualitatively reproduce the inhibiting effect on tumor growth caused by the tumor mass effect and a history of BPH before PCa detection. To verify the validity of the use of linear kinematic theory, we computed the symmetrical and skew-symmetrical components of the displacement gradient with respect to the initial configuration and determined that they were sufficiently small that their products were not substantial compared with the linear terms and, in particular, that the skew-symmetric components were negligible everywhere. Consequently, we feel confident that the computed results provide physically meaningful information.

Numerical Methods.

Since tumor growth and BPH are modeled as quasistatic processes, we adopted a staggered approach to solve the equations in our mathematical model. We calculated tumor growth at every time step, but we updated the displacements only every two time steps. We performed spatial discretization by means of a standard isogeometric Bubnov–Galerkin approach using a 3D C1 quadratic nonuniform rational B-spline (NURBS) space (39, 40, 60⇓–62). Temporal integration in the tumor growth problem was carried out with the generalized-α method (63, 64). This technique led to a nonlinear problem in each time step, which we linearized using the Newton–Raphson method. The resulting linear system was solved using the generalized minimal residual method (GMRES) algorithm (65) with a diagonal preconditioner. We also used the preconditioned GMRES algorithm to solve the quasistatic elastic problem. We chose a constant time step of 0.002 y and the parameters in the generalized-α method were set as in previous studies (40, 66).

Visualization.

We used ParaView (67) to visualize and explore the results of our simulations. We used the isovolume ϕ≥0.5 to represent the tumor, which permitted us to study the evolution of PCa morphologies. From the C1-continuous basis, we computed the stress fields pointwise. We identified regions of interest for the stress fields in the simulations and isolated them using appropriate geometric filters in ParaView. The stresses reported in the text correspond to the range of values that best describe the stress field in each region of interest.

Construction of the Prostate Mesh.

Multiple methodologies permit the construction of solid anatomic NURBS models (68, 69). Because the geometries of the prostate and a solid torus are topologically equivalent, we leveraged a parametric mapping algorithm (70, 71) to deform a torus solid NURBS model to match with a patient-specific prostate surface model. We used 3DSlicer (72) to generate a triangular surface model of the prostate from the contours of the organ and the urethra drawn on the T2-weighted MR images, using the provided prostate segmentation as guidance. The resulting surface was smoothed in MeshLab (73).

The original torus and prostate NURBS meshes were discretized with 32 × 32 × 8 elements along the toroidal direction, the cross-section circumferential direction, and the cross-section radial direction, respectively. We globally refined the prostate mesh to 256 × 256 × 64 elements, using standard knot insertion (61) to perform our simulations with a good level of accuracy.