A Pilot Study on Patient-specific Computational Forecasting of Prostate Cancer Growth during Active Surveillance Using an Imaging-informed Biomechanistic Model

Prostate cancer is the second most common type of cancer and the fifth leading cause of cancer death among men worldwide (1). The clinical management of prostate cancer relies on two key strategies (2): regular screening and patient triage into risk groups. Regular screening of prostate cancer is performed in men over 50 and consists of a digital rectal exam and the measurement of the serum level of the PSA. If the results of these two tests are suggestive of prostate cancer, then a multiparametric MRI (mpMRI) scan will usually be performed to assess the prostate noninvasively and also guide a subsequent biopsy for histopathologic confirmation (2, 3). On the basis of the data collected during this diagnostic stage, clinical protocols assign risk groups and associate specific management options to each of them (2). The definition of these risk groups usually relies on the serum PSA, the clinical stage (i.e., an estimate of tumor size and extent) according to the tumor-node-metastases scale, and the biopsy Gleason score (GS), which is the gold standard histopathologic marker of prostate cancer aggressiveness and is defined upon the cancerous tissue architecture observed in biopsy samples or surgically-removed prostates (2, 4). Thanks to the current clinical management of prostate cancer, most tumors are detected and successfully treated at early organ-confined stage, which usually poses a low to intermediate risk to the patient (5, 6). Although standard treatments for prostate cancer (e.g., surgery, radiotherapy) exhibit high rates of overall and disease-free survival (2, 7, 8), many newly diagnosed prostate cancer cases are known to be indolent and may not produce any symptoms or require treatment for a long time. These patients are prone to potential treatment side effects (e.g., incontinence, impotence) that can adversely impact their quality of life without improving longevity. Thus, overtreatment of prostate cancer is a major concern in the current clinical management of the disease (3, 9, 10). In addition, undertreatment of prostate cancer due to limitations in diagnosis and management (e.g., biopsy sampling errors, mpMRI confounders, limited age and comorbidity assessment) constitutes another important clinical challenge, because it may result in rapid growth of aggressive tumors, treatment failure, and reduced survival (3, 9–11).

As an alternative to direct treatment after diagnosis, active surveillance (AS) is another standard-of-care management option that is suitable for many newly diagnosed prostate cancer cases exhibiting low to intermediate clinical risk (2, 3, 8, 12). In AS, patients are closely monitored via longitudinal mpMRI scans, digital rectal exams, PSA tests, and biopsies, such that treatment is delayed until these monitoring tests confirm prostate cancer progression to higher risk. Thus, AS is widely regarded as an ideal clinical strategy to reduce the overtreatment of newly diagnosed prostate cancer (2, 3, 12). In addition, the changes in longitudinal data collected during AS (e.g., mpMRI-observable features, PSA kinetics) can be leveraged to inform clinical decision-making and, hence, contribute to reduce treatment excesses and deficiencies (12, 13). However, current AS protocols largely rely on recommending monitoring tests either at fixed times (e.g., mpMRI scan potentially followed by biopsy every 6–36 months, PSA test every 3–12 months) or after the observation of clinical events suggesting tumor progression (e.g., imaging after a sustained PSA rise, biopsy after radiological worsening of an mpMRI-detected lesion; see refs. 2, 12, 13). This population-based, observational approach largely ignores the complex underlying tumor dynamics and the intrinsic heterogeneity of prostate cancer between and within patients (3, 4, 13–15), complicates the design of personalized plans to ensure optimally controlled tumor monitoring, and hinders the early detection of prostate cancer progression, which may occur between two consecutive and widely-spaced monitoring tests. Thus, the optimal timing and type of monitoring tests for each individual prostate cancer case constitutes a central challenge in AS, along with the definition of triggering criteria to switch from monitoring to radical treatment (2, 3, 12, 13, 16). To address these unmet and timely demands in AS, we propose to leverage personalized computational predictions of prostate cancer growth (17, 18).

Computational tumor forecasting is a technology that enables the estimation of the growth of tumors and their response to treatment on a personalized basis (18). These tumor forecasts are obtained by leveraging mathematical models based on differential equations or agent-based formulations that describe the main biological mechanisms underlying the development of these diseases and the effect of specific treatments (18, 19). By calibrating the model parameters governing tumor dynamics with patient-specific data, the ensuing personalized model enables the representation of observed tumor growth as well as the prediction of the future development of the tumor, which can be leveraged to guide monitoring strategies or adjust treatment regimens for each individual patient (18–26). Indeed, the increasing success of computational tumor forecasting has motivated its use to design digital twin technologies to systematically guide clinical decision-making using a predictive, patient-specific approach (27). In particular, medical imaging data, such as mpMRI scans, provide spatially-resolved anatomic and physiologic information about tumors that enable the personalization of spatiotemporal biomechanistic models based on partial differential equations, which can then render three-dimensional (3D) tumor forecasts within the patient's affected organ (18). While this type of models has been extensively used to predict the development and therapeutic response of other tumors (20–23, 28–31), there is a dearth of applications for prostate cancer (17, 32–35). Instead, most biomechanistic models developed for prostate cancer consist of time-resolved differential formulations primarily based on serum PSA dynamics that have been successfully used to predict response to specific treatments, especially in the advanced stage of the disease (24–26, 36). The abundance of these PSA models is primarily motivated by the extensive use of this biomarker to monitor prostate cancer, while longitudinal mpMRI data are almost only available in AS. However, time-resolved models lack spatial information on prostate cancer growth that can be key during AS, for example, to plan biopsies or to assess malignancy in prostate areas that may associate with worse prognosis (e.g., prostate capsule, regions adjacent to neurovascular bundles and seminal vesicles; see refs. 2–4, 13).

Here, we present a pilot study of personalized forecasting of untreated prostate cancer growth during AS using a spatiotemporal biomechanistic model informed by longitudinal mpMRI data collected for each individual patient. Toward this end, we develop a clinical-computational pipeline to segment the prostate and tumor, perform intrascan and interscan registration of imaging data, build a 3D virtual representation of the prostate geometry holding the longitudinal spatial changes of the tumor, calibrate a biomechanistic model of prostate cancer growth, and perform patient-specific tumor forecasting (see Fig. 1). First, we analyze whether our computational framework can represent prostate cancer growth by informing the model with three imaging datasets for each patient. Then, we investigate whether we can predict prostate cancer growth after informing the model with two imaging datasets by comparing a prediction of the personalized model against a third mpMRI dataset. In addition, because prostate cancer risk and AS eligibility heavily rely on GS (2, 12), we further analyze our personalized computer simulations to search for potential model-based predictive biomarkers of prostate cancer progression to higher-risk disease.

Anonymized patient data were retrospectively collected at The University of Texas Health Science Center San Antonio (San Antonio, TX) following an Institutional Review Board (IRB)-approved and Health Insurance Portability and Accountability Act–compliant protocol, which did not require informed consent. The inclusion criteria were: (i) newly diagnosed, untreated, organ-confined prostate cancer managed with AS after diagnosis; (ii) availability of three mpMRI scans during AS; (iii) a visible tumor across all mpMRI scans; and (iv) histopathologically-confirmed prostate cancer via biopsy and/or analysis of the surgically-removed prostate after AS. The cohort leveraged in this study consisted of 7 patients (n = 7), each of them having one organ-confined tumor matching the inclusion criteria. The longitudinal mpMRI scans of each patient were collected over a period of median (range) of 3.4 (1.4, 4.9) years with interscan frequency of 1.4 (0.5, 2.6) years. For each patient, we define a relative time scale beginning at the time of the first mpMRI scan (t = 0). Although most biopsies were performed after an mpMRI scan, some patients also had biopsies before the first imaging session. Baseline GS at t = 0 was assigned according to the GS of the first biopsy performed within 3 months after the first mpMRI date. If there was no biopsy within this period but the patient had a posterior biopsy (i.e., >3 months after the first mpMRI scan) and a preimaging biopsy, the first postimaging GS was used as the baseline measurement only if it was identical to the preimaging GS value. If the patient had no biopsies after the first mpMRI scan, the most recent GS measured before the first mpMRI was used as baseline. All other postimaging GS measurements were used for temporally-informed clinical risk stratification during AS. This process enabled the definition of GS at t = 0 for all patients and resulted in a total of 16 eligible GS measurements with a median (range) of 3 (1, 4) values per patient. These GS measurements were evenly distributed in two subgroups exhibiting a GS = 3 + 3 (n = 8) and a GS ≥ 3 + 4 (n = 8), which we denote as lower-risk and higher-risk prostate cancer. This division enabled a balanced subgroup definition and is consistent with the fact that prostate cancer with GS 3 + 3 is extensively recognized as eligible for AS, while there is debate regarding the inclusion in AS protocols of prostate cancer cases exhibiting Gleason grade 4 (2, 12).

The mpMRI datasets were collected using 3T scanners for all subjects. A TIM Trio scanner (Siemens Healthineers USA) housed at the Research Imaging Institute of The University of Texas Health Science Center San Antonio (San Antonio, TX) was used to obtain a baseline scan with a cardiac coil used on the pelvis (37, 38). Follow-up scans were performed following the standard-of-care using either a GE 750W (GE Healthcare USA) at the South Texas Veterans Healthcare System or a Siemens Skyra scanner (Siemens Healthineers USA) at The University of Texas Health Science Center San Antonio (San Antonio, TX). The mpMRI acquisition protocols included T₁-weighted, T₂-weighted, diffusion-weighted, and dynamic contrast-enhanced sequences (T1W, T2W, DW, and DCE, respectively). The main acquisition parameters were as follows: echo time <90 seconds, repetition time >3,000 ms, slice thickness <4.0 mm without gap, field of view of 160–220 mm, and b-values of approximately 50, 400, 800, and 1,400 s/mm² on all scanners. Apparent diffusion coefficient (ADC) maps were calculated from DW-MRI data using the standard monoexponential model (39).

All biopsies were performed using a transrectal ultrasound approach. After an mpMRI was performed, the regions of interest and prostate boundaries were marked in DynCAD software as part of the UroNav (Phillips Healthcare) biopsy system. During biopsy, once the prostate was scanned and registered with the MRI, three targeted biopsies were obtained from each region of interest and a standard 12-core prostate biopsy was also completed. Biopsy specimens underwent standard pathologic processing and were graded using the International Society of Urologic Pathology criteria and standard Gleason grade scoring. The same histopathologic assessment process was used for standard biopsies before the first mpMRI scan and for surgically-removed prostate specimens.

The integration of the mpMRI data from each patient within our biomechanistic model of prostate cancer growth required five preprocessing operations: (i) segmentation of the prostate in each T2W image, (ii) intrascan registration of the T2W map to the ADC map for each mpMRI scan, (iii) segmentation of the tumor on each ADC map, (iv) interscan registration of the longitudinal mpMRI data and the segmentations for each patient, and (v) conversion of intratumoral ADC values to tumor cell density estimates. This computational pipeline was designed following previous efforts developed to support mpMRI-informed personalized forecasting for brain (29) and breast cancers (21, 28). For each mpMRI dataset, the segmentation of the prostate and the tumor, as well as the intrascan registration, were carried out using a combination of manual and semi-automatic tools available in 3DSlicer (40). Toward this end, the T2W and ADC maps were first resampled to isotropic resolution using linear interpolation. We next briefly describe the preprocessing operations in the order they were implemented.

The prostate gland was manually segmented in the resampled T2W images using anatomic landmarks as reference (3, 4). To facilitate the ensuing intrascan registration, we also segmented the prostate in the ADC maps.

To align the ADC map to the T2W map of each mpMRI dataset, we performed a rigid intrascan registration by leveraging the General Registration module of 3DSlicer and using the corresponding prostate segmentations for masking. Following this intrascan registration, we smoothed the T2W-based segmentations of the prostate using the Segmentation smoothing module in 3DSlicer.

The segmentation of the tumor consisted of a multistep, hierarchical procedure. First, we manually drew a gross segmentation of the tumor over the registered ADC map based on the tumor location provided in the mpMRI reports, which was confirmed by an experienced urologist who specializes in prostate cancer and its detection using mpMRI (M.A. Liss). Second, we manually delineated a volume of healthy tissue in the same region of the prostate where the tumor is primarily located (i.e., peripheral zone or central gland) with approximately the same volume as the gross segmentation of the tumor. We then calculated the mean ADC over this volume as an estimate of ADC in healthy tissue (ADC_h) using the Segment Statistics module in 3DSlicer. Third, we used the threshold tool in 3DSlicer to delineate the tumor core within the gross segmentation of the tumor, which was defined as the isovolume of ADC values with less than 70% of ADC_h. This value was selected after an analysis of five studies investigating the correlation of ADC and GS (41–45). These studies show that prostate cancer exhibits lower ADC than healthy prostatic tissue and that ADC progressively decreases as the GS increases. This trend seems to plateau for higher GS values, which aligns with a higher tumor cell density approaching tissue carrying capacity in those tumors (46–48). By fitting a hyperbolic tangent function to the average of the ADC ratios (ADC/ADC_h) obtained from the mean ADC_h and mean GS-specific ADC values in the five studies (41–45), we determined that the 70% of ADC_h was a sufficient threshold to identify mpMRI-visible tumors (see Fig. 2A). This fitting was performed using the Curve Fitting Toolbox in MATLAB (R2021b; The Mathworks) and defining a continuous GS scale ranging from 0 to 10, where values between 0 and 2 represent healthy tissue and pretumoral lesions (4). In addition, the choice of the hyperbolic tangent includes a slow decrease of ADC in pretumoral lesions and low GS tumors, which makes them practically indistinguishable from healthy tissue and aligns with the fact that mpMRI-observable prostate cancer generally exhibits GS ≥ 3 + 3 at diagnosis (3, 14, 16). Supplementary Data S1 provides further details about this hyperbolic tangent fit. Finally, the resulting tumor core was expanded to define the final tumor segmentation by leveraging the margin tool in 3DSlicer. This technique was applied within the gross tumor segmentation and using a margin of 2 to 4 mm (equivalent to a kernel from 3 × 3 × 3 to 5 × 5 × 5 voxels, depending on the tumor size), such that the tumor border lied over larger ADC values (e.g., ≳90% of ADC_h). Hence, this margin facilitated a smooth transition of ADC values from the tumor core to the neighboring healthy tissue that matches the natural solution of our biomechanistic model.

For each patient, the three mpMRI datasets collected during AS and their corresponding segmentations were co-registered to a common data frame, which was set as that of the first mpMRI scan (see Fig. 1C). This interscan registration was carried out using a biomechanically-constrained, deformable image registration algorithm to control for bulk deformations of the prostate anatomy between scan dates while preserving imaging features associated with longitudinal disease-related changes. The registration process consisted of a rigid registration of the T2W-based prostate segmentations followed by a deformable linearized iterative boundary reconstruction algorithm (49, 50). This registration method filters out purely elastic deformation effects between imaging dates while preserving the effect of nonelastic imaging changes associated with tumor growth during AS. Further details are available in Supplementary Data S2. The interscan registration was applied to the T2W data, ADC maps, and tumor segmentations of each patient. After the interscan registration, we updated the reference prostate segmentation (i.e., that defined in the T2W image of the first mpMRI scan) to prepare it for integration within our modeling framework by subtracting the urethral region. To do this, we used the eraser tool in 3DSlicer with a spherical geometry and 4-mm diameter. Anatomic landmarks in the T2W data were used for reference.

Because ADC has been shown to decrease as tumor cell density increases in higher GS tumors (41–48), the ADC values within each tumor segmentation were converted to normalized tumor cell density by leveraging a formulation that has been successfully employed to represent the inverse relationship between ADC and tumor cell density in other solid tumors (18, 21, 28, 29), as follows:

In Eq. (1), (x, t) is the normalized tumor cell density at position x and time t (i.e., 0 ≤ (x, t) ≤ 1), which is defined as the ratio of tumor cell density N(x, t) to the tissue carrying capacity θ (i.e., the maximally admissible tumor cell density in the prostate, 0 ≤ N(x, t) ≤ θ). In addition, ADC_h is the ADC in the healthy tissue as defined above, ADC(x, t) is the ADC at position x and time t within the tumor, and ADC_min is the minimum ADC observable within a tumor. For each ADC map, the latter was estimated from the lower horizontal asymptote of the hyperbolic tangent function in Fig. 2A as ADC_min/ADC_h = 0.25. Figure 2B further provides the mapping through Eq. (1) of the hyperbolic tangent function describing the changes in ADC ratio with respect to GS, thereby yielding the changes in normalized tumor cell density with respect to GS within our modeling framework. Because ADC_h and ADC_min are fixed estimates, ADC values above ADC_h and below ADC_min were truncated to render a normalized tumor cell density of 0 and 1, respectively. This operation avoided unphysical negative tumor cell densities and values over the carrying capacity. Thus, according to Eq. (1) and the plots in Fig. 2, the segmentation of the tumor core using values below 70% ADC_h corresponds to normalized tumor cell density values ≥ 0.4, and, after the margin expansion, the border of the resulting tumor segmentation has values ≲ 0.15.

We model prostate cancer growth during AS in terms of the spatiotemporal dynamics of the tumor cell density, N = N(x, t), using a reaction-diffusion partial differential equation, as follows:

Equation (2) is commonly known as the Fisher–Kolmogorov equation (51). Many patient-specific tumor forecasting studies have directly relied on this equation or leveraged it as a basis to build more sophisticated models (18, 21–23, 28–31). The right-hand side of Eq. (2) describes the dynamics of the tumor cell density as a combination of two mechanisms: tumor cell mobility, which is represented by a diffusion process governed by the tumor cell diffusivity (D), and tumor cell net proliferation, which is modeled with a logistic term controlled by the net proliferation rate (ρ) and the tissue carrying capacity (θ). The tumor cell net proliferation rate encompasses the balance of tumor cell proliferation and death. In this work, we assume that the carrying capacity is constant for each patient, which enables posing our model in terms of the normalized tumor cell density ( = N/θ) by dividing both sides of Eq. (2) by θ:

Hence, once we calibrate the values of parameters D and ρ in Eq. (3) that best explain prostate cancer dynamics as observed in the longitudinal mpMRI measurements of normalized tumor cell density collected during AS for each patient, we can obtain personalized forecasts of tumor growth. The rationale for adopting the model formulation in Eq. (3) is that we want to leverage a parsimonious formulation of prostate cancer growth, such that it features a minimum set of personalized parameters that does not require extensive datasets, facilitates parameter identifiability, and avoids overfitting. Although θ may vary for each patient in this modeling framework, the normalization underlying Eq. (3) enables us to focus on relative changes of tumor cell density that match the changes in ADC observed in prostate cancer with increasing GS values in the literature (41–45). Thus, in our modeling framework θ can be patient-specific, but we do not need it to interrogate spatiotemporal prostate cancer dynamics. Figure 3 illustrates how tumor cell mobility and net proliferation contribute to prostate cancer dynamics according to the models in Eqs. (2) and (3). As shown in Fig. 3, these models are posed within the patient's prostate geometry, which corresponds to the T2W-based segmentation from the first mpMRI scan. Furthermore, the model is defined over a temporal interval [0, T], where T is a specific time horizon (e.g., for model calibration or tumor forecasting for each patient). To complete the model, we need to define boundary conditions over the prostate surface. Because the prostate cancer cases eligible for this study were organ-confined, we set zero-flux boundary conditions (i.e., ∇N ∙ n = 0 or, equivalently, ∇ ∙ n = 0; where n is an outward unit vector orthogonal to the prostate boundary). In addition, the initial conditions for Eq. (3), (x, 0), were always defined with the normalized tumor cell density map obtained from the first mpMRI scan.

To calibrate the model and obtain prostate cancer predictions for each patient, we performed computer simulations of our model by solving Eq. (3) numerically. These computer simulations constitute an in silico representation of personalized prostate cancer growth over the patient's prostate anatomy and the temporal interval [0, T], which we can compare with mpMRI data at scan times. We leveraged isogeometric analysis (IGA), which is a high-fidelity generalization of the finite element method that uses highly smooth functions to discretize the geometry and the solution of partial differential equations (18, 52). In particular, we discretized Eq. (3) in space using a standard isogeometric Bubnov–Galerkin method that relied on a 3D C¹ quadratic non-uniform rational B-spline (NURBS) functional space (18, 33, 52). This method requires the construction of a patient-specific mesh of the prostate using the same function space. Toward this end, we used a parametric mapping method (17, 32), which deformed a reference torus mesh to match the final prostate segmentation defined over the T2W image from the first mpMRI scan of each patient. The original torus and prostate meshes had a polar discretization with 32 elements in the circumferential directions and eight elements in the radial direction. Each prostate mesh was subsequently refined using standard knot insertion (52) for numerical accuracy in the computer simulations of our biomechanistic model. Hence, the final patient-specific prostate meshes had 128 elements in the circumferential directions and 32 elements in the radial direction. The tumor segmentations and corresponding normalized tumor cell density maps were then L²-projected over the prostate mesh for each patient (17, 32, 33). Furthermore, the time interval [0, T] was discretized with a constant time step Δt = 1 day for all patients, which sufficed to capture the spatiotemporal dynamics observed in the patients of this study. We integrated in time using the generalized-α method (33, 52). This iterative algorithm produces a series of nonlinear systems of equations in each time step, which we linearized with the Newton-Raphson method. The ensuing linear systems were solved utilizing the generalized minimal residual method. Further details on the implementation of these numerical methods to solve our biomechanistic model are available in refs. 17, 18, 32, 33, 52.

The personalized parameterization of our biomechanistic model consisted of finding the values of the parameter set (D, ρ) that minimize the relative squared difference between the model predictions and the corresponding imaging-based measurements of normalized tumor cell density at the times of the mpMRI scans used to inform the model. This calibration problem can be formulated as:

where n_s is the number of mpMRI scans used for model calibration, Ω represents the 3D geometry of the patient's prostate, and and denote the normalized tumor cell density measured from imaging data and calculated from the model with parameters D and ρ at the scan times t_i (i = 1, …, n_s), respectively. Hence, Eq. (4) defines a nonlinear least-squares problem to find the patient-specific parameter values (D, ρ), which we solved by leveraging the Gauss–Newton method (53). The initial guess for the parameters was set as D = 5 · 10⁻³ mm²/day and ρ = 2 · 10⁻³ 1/day, while the admissible values of the parameters were constrained to D ∈ [1.0 · 10⁻⁶, 10] mm²/day and ρ ∈ [1.0 · 10⁻⁶, 1] 1/day.

Model performance during calibration and forecasting was assessed with a panel of global and local metrics. The global metrics were the tumor volume (V_T) and the total tumor cell volume (V_N), which are defined as:

In Eqs. (5) and (6), Ω_T is the 3D tumor region where with being a threshold value to systematically identify tumor tissue both in the imaging data and the model forecasts. Hence, this threshold facilitated the automatic delineation of a tumor region in the imaging data, for which the values over the tumor segmentation border were ≲ 0.15, and from the model simulations, which provided a continuous map within the range [0,1]. For this purpose, we set _th = 0.15 in all the calculations presented in this work. The global metrics V_T and V_N were calculated patientwise in each calibration scenario, and model-data agreement was assessed via the Pearson and concordance correlation coefficients over the cohort (PCC and CCC, respectively). The local metrics employed to analyze the performance of the model were: the Dice similarity coefficient (DSC) as well as the root mean squared error (RMSE), the local PCC, and the local CCC between the tumor cell density maps calculated from the imaging data and leveraging our biomechanistic model. These local metrics were calculated patientwise and we analyzed their distribution using boxplots.

We analyzed the personalized computer simulations of our model to find model-based biomarkers of high-risk prostate cancer. Toward this end, we first calculated the values of a panel of six candidate markers at the times of histopathologic assessment of the tumors in our patient cohort, and then compared their values between the subgroups of low-risk and high-risk prostate cancer (i.e., GS = 3 + 3 and GS ≥ 3 + 4, respectively). These candidate markers were calculated using the model simulations obtained with the personalized parameters from a global calibration study, in which the model was informed with the maximum amount of mpMRI data available for each patient in this study. The list of candidate markers was designed on the basis of previous tumor forecasting studies (22, 23, 28–32, 54) and our analysis of the personalized predictions in this work. The six candidate markers were: the prostate volume (V_P), the tumor volume (V_T), the total tumor cell volume (V_N), the mean normalized tumor cell density (), the total tumor index (N_T), and the mean proliferation activity of the tumor (A_p). The prostate volume and the last three markers were calculated as follows:

The prostate volume and the total tumor index defined in Eqs. (7) and (9), respectively, were included to account for the organ size effect on tumor dynamics, because high-volume prostates enlarged by concomitant benign prostatic hyperplasia (BPH) have been suggested to exhibit a lower risk of high-grade prostate cancer (32, 55). In addition, the mean proliferation activity of the tumor (A_p) was defined on the basis of the work presented in ref. 54.

Any candidate marker exhibiting a statistically significant difference between the low-risk and high-risk prostate cancer subgroups was considered an eligible model-based biomarker of high-risk prostate cancer. We analyzed the performance of these biomarkers by building univariate logistic regression classifiers and using ROC curve analysis. We further built bivariate logistic regression classifiers including a maximum of two markers from the candidate list, such that at least one was identified as a biomarker. We did not consider logistic regression classifiers with more than two input markers to avoid overfitting given the reduced size of the histopathologic measurements in the cohort (n = 16). All logistic regression classifiers were built using the function glmfit from the Statistics and Machine Learning Toolbox in MATLAB. The area under the ROC curve was calculated numerically with the trapezoidal rule and the optimal performance point was determined using the minimum distance to the upper left corner (i.e., where sensitivity and specificity are maximal).

Given the small size of the cohort in this work (n = 7), we did not partition it to define a patient subgroup to train the classifiers and another group to validate them. Nevertheless, to investigate their predictive performance, we analyzed them using the personalized computer simulations from the fitting-forecasting study. In particular, we analyzed whether the best performing classifier, which was trained with markers calculated with personalized simulations from the global calibration study, was able to anticipate the development of high-risk disease for each patient when the model was only informed by the first two mpMRI scans.

We leveraged Wilcoxon rank-sum and signed-rank tests to compare the local metrics of model performance calculated at the times of the second and third mpMRI scans in both the global calibration scenario and the fitting-forecasting study. These two statistical tests were also used to investigate differences in the model parameters as well as global and local metrics of model performance between the global calibration scenario and the fitting-forecasting study. In addition, the Wilcoxon rank-sum test was also used to identify differences between the candidate markers of high-risk prostate cancer between the low-risk and high-risk prostate cancer subgroups. The two types of statistical tests were performed using the functions ranksum and signrank from the Statistics and Machine Learning Toolbox in MATLAB. In the Results section, we specify when we use each type of test and whether it is two-tailed or one-tailed for each statistical analysis. The level of significance for all statistical tests was set to 5%.

The raw data leveraged in this study were generated at The University of Texas Health Science Center San Antonio (San Antonio, TX). Raw data are not publicly available under the IRB data usage agreement for this study, but the derived data that support the findings presented in this work are available from the corresponding author upon reasonable request.

We first performed the global calibration study to assess the ability of our computational framework to represent untreated prostate cancer growth during AS. In this study, the biomechanistic model was informed by the three mpMRI scans available for each patient (n = 7). The reference prostate geometry from the first mpMRI scan had median (range) volume of 36.4 (18.5, 67.3) cc. The imaging-based segmentation of the tumors over the first, second, and third mpMRI datasets resulted in a tumor volume (V_T) of 0.12 (4.4 · 10⁻³, 1.26) cc, 0.28 (0.06, 4.44) cc, and 0.49 (0.14, 7.12) cc, respectively. The corresponding total tumor cell volumes (V_N) were 0.07 (1.5 · 10⁻³, 0.91) cc, 0.13 (0.02, 2.45) cc, and 0.31 (0.06, 3.89) cc. The global calibration study resulted in a distribution of tumor cell diffusivity coefficients (D) and net tumor cell proliferation rates (ρ) of 1.26 · 10⁻³ (5.74 · 10⁻⁴, 5.35 · 10⁻³) mm²/day and 1.86 · 10⁻³ (1.09 · 10⁻³, 5.13 · 10⁻³) 1/day, respectively. After initializing our biomechanistic model with the normalized tumor cell density maps from the first mpMRI scan and corresponding calibrated parameters for each patient, the ensuing personalized simulations resulted in tumor volumes (V_T) of 0.33 (0.04, 5.44) cc and 0.59 (0.15, 7.69) cc at the time of the second and third mpMRI scans, respectively. The corresponding model-calculated total tumor cell volumes (V_N) were 0.13 (0.01, 2.28) cc and 0.23 (0.04, 3.37) cc.

Figure 4 illustrates the 3D imaging measurements and calculations of the tumor geometry and the normalized tumor cell density map obtained with the personalized model for 4 patients at the second and third imaging timepoints. The corresponding results for the other 3 patients are provided in Supplementary Fig. S1. In addition, Fig. 5A and B provide unity plots comparing the imaging measurements and the model calculations of the tumor volume (V_T) and the total tumor cell volume (V_N). Considering the second mpMRI scan for each patient, the PCC and CCC for V_T were 0.99 and 0.97, respectively. The corresponding PCC and CCC for V_N were 0.99 and 0.98. At the time of the third mpMRI scan, the PCCs and CCCs of both V_T and V_N were ≥0.99 in all cases. Figure 5C–F further provide the distribution of the local metrics assessing the pointwise agreement between the imaging-based and the model-calculated normalized tumor cell density maps ((x, t)) across the patient cohort (n = 7). At the second imaging timepoint, we obtained a DSC of 0.81 (0.66, 0.85), an RMSE of 0.022 (0.008, 0.071), a local PCC of 0.50 (0.19, 0.68), and a local CCC of 0.49 (0.18, 0.67). The corresponding distributions of DSC, RMSE, local PCC, and local CCC at the time of the third mpMRI scan were 0.77 (0.64, 0.85), 0.044 (0.016, 0.107), 0.28 (0.18, 0.63), and 0.28 (0.17, 0.61). No significant differences were detected between the values of each local metric at the second and third imaging timepoints under Wilcoxon rank-sum testing (P > 0.05). Nevertheless, two-sided Wilcoxon signed-rank tests at the time of the third mpMRI scan detected significantly higher RMSE (P = 0.016) and significantly lower DSC, local PCC, and local CCC (P = 0.047 in all three cases).

To investigate the predictive performance of our computational framework, we carried out a fitting-forecasting study: for each patient, we first initialized our biomechanistic model with the data from the first mpMRI, we then calibrated the parameters using the data from the second mpMRI dataset, and we finally validated a personalized forecast of prostate cancer growth against the imaging data collected at the third imaging timepoint. The median (range) of the tumor cell diffusivity coefficients (D) and net tumor cell proliferation rates (ρ) obtained in the fitting-forecasting study for each patient (n = 7) were 1.40 · 10⁻³ (5.35 · 10⁻⁴, 4.18 · 10⁻³) mm²/day and 2.41 · 10⁻³ (7.84 · 10⁻⁴, 5.15 · 10⁻³) 1/day, respectively. The ensuing patient-specific model simulations resulted in tumor volumes (V_T) of 0.29 (0.04, 4.98) cc and 0.45 (0.16, 6.95) cc at the second and third imaging timepoints, respectively. The corresponding total tumor cell volumes (V_N) were 0.12 (0.01, 2.20) cc and 0.20 (0.05, 3.52) cc.

Figure 6 presents the model-predicted tumor geometry and normalized tumor cell density map obtained in the fitting-forecasting study at the times of the second and third mpMRI scans for the same 4 patients shown in Fig. 4 in the global calibration study, along with the corresponding imaging-based measurements. Similar results for the other 3 patients are provided in Supplementary Fig. S2. In addition, the unity plots in Fig. 7A and B compare the model predictions of the tumor volume (V_T) and the total tumor cell volume (V_N) to their corresponding imaging measurements. At the time of the second mpMRI scan (i.e., the calibration horizon), the PCC and CCC values for the model-data agreement of V_T and V_N were ≥0.99 in all cases. The PCC and the CCC for V_T resulted in 0.97 and 0.96 at the third imaging timepoint (i.e., the forecasting horizon), respectively. The corresponding values of the PCC and the CCC for V_N at the third mpMRI date were 0.93 in both cases. Figure 7C–E also present the distributions of the local metrics quantifying the pointwise agreement between the imaging measurements and model forecasts of the normalized tumor cell density maps ((x, t)) over the patient cohort (n = 7). At the second imaging timepoint, we obtained a DSC of 0.82 (0.67, 0.85), an RMSE of 0.022 (0.008, 0.068), a local PCC of 0.55 (0.19, 0.67), and a local CCC of 0.55 (0.18, 0.66). The corresponding values of the DSC, RMSE, local PCC, and local CCC at the forecasting horizon were 0.76 (0.65, 0.85), 0.048 (0.016, 0.137), 0.37 (0.26, 0.61), and 0.33 (0.24, 0.60), respectively. Two-sided Wilcoxon rank-sum testing resulted in no significant differences (P > 0.05) between the values of the local metrics at the times of the second and third mpMRI scans, although the DSC values were identified as significantly higher at the calibration horizon with respect to the forecasting timepoint under a one-sided Wilcoxon rank-sum test (P = 0.049). According to two-sided Wilcoxon signed-rank tests, model forecasts at the third imaging timepoint exhibited a significantly higher RMSE (P = 0.016) and significantly lower DSC, local PCC, and local CCC (P = 0.016, 0.047, and 0.047, respectively).

No significant differences were detected between the values of the model parameters (D, ρ) obtained in the global calibration and the fitting-forecasting study under both Wilcoxon rank-sum and signed-rank testing (P > 0.05). Nevertheless, we observed a tendency toward lower tumor cell diffusivity coefficients (D) and higher net tumor cell proliferation rates (ρ) in the fitting-forecasting study for each patient under one-sided Wilcoxon signed-rank tests (P = 0.055). In addition, no Wilcoxon rank-sum test identified significant differences between the local model-data agreement metrics from the global calibration and the fitting-forecasting studies at either the second or third imaging timepoints (P > 0.05). However, two-sided Wilcoxon signed-rank testing resulted in significantly lower RMSE values at the third mpMRI date in the global calibration study (P = 0.016) and at the second imaging timepoint in the fitting-forecasting study (P = 0.016). Furthermore, the DSC at the second scan timepoint was significantly higher in the fitting-forecasting study under a one-sided Wilcoxon signed-rank test (P = 0.039).

To search for potential model-based biomarkers of high-risk prostate cancer, we investigated six quantities of interest that were calculated at the times of histopathologic assessment of each patient's tumor (n = 16) using the personalized model simulations from the global calibration study. The median (range) of the prostate volume (V_P), the tumor volume (V_T), the total tumor cell volume (V_N), the mean normalized tumor cell density (), the total tumor index (N_T), and the mean proliferation activity of the tumor (A_p) in the low-risk prostate cancer group (i.e., GS = 3 + 3; n = 8) were 37.7 (21.8, 39.5) cc, 0.17 (0.01, 5.09) cc, 0.07 (1.5 · 10⁻³, 2.12) cc, 0.39 (0.26, 0.52), 2.31 · 10⁻³ (3.82 · 10⁻⁵, 5.36 · 10⁻²), 3.61 · 10⁻⁴ (2.28 · 10⁻⁴, 4.94 · 10⁻⁴) 1/day. The corresponding median (range) values in the high-risk prostate cancer subgroup (i.e., GS ≥ 3 + 4, n = 8) were 34.0 (18.5, 67.3) cc, 0.58 (0.12, 5.90) cc, 0.19 (0.05, 2.58) cc, 0.39 (0.30, 0.54), 9.61 · 10⁻³ (1.69 · 10⁻³, 6.53 · 10⁻²), 3.88 · 10⁻⁴ (3.57 · 10⁻⁴, 1.07 · 10⁻³) 1/day. Figure 8A–F plot the distribution of these candidate markers in the low-risk and high-risk prostate cancer subgroups. Only the mean proliferation activity of the tumor (A_p) was significantly larger in high-risk prostate cancer according to a one-sided Wilcoxon rank-sum test (P = 0.041). Figure 8G shows the ROC curve for a univariate logistic classifier of high-risk prostate cancer that was constructed using the mean proliferation activity of the tumor (A_p). The AUC of this ROC curve is 0.77 and the optimal performance point operates at 75% sensitivity and specificity. In addition, Fig. 8G also shows the ROC curve for the best performing bivariate logistic classifier that could be built using the mean proliferation activity of the tumor (A_p) and one of the other markers, which resulted to be the total tumor index (N_T). The AUC for this ROC curve is 0.83 and the optimal performance point also operates at 75% sensitivity and specificity.

To investigate whether the bivariate logistic classifier would anticipate the detection of high-risk prostate cancer at the time of the second mpMRI scan (i.e., once the biomechanistic model is fully calibrated for each patient, and we can perform a personalized forecast), we tested it in the fitting-forecasting scenario. The results of this analysis for the 4 patients in Figs. 4 and 6 are shown in Fig. 9, while the corresponding results for the other 3 patients are reported in Supplementary Fig. S3. These figures further compare the performance of the classifier using the personalized model simulations from the global calibration scenario (i.e., those used for classifier training). The prostate cancer cases in Fig. 9A and D were consistently classified as low-risk and high-risk prostate cancer, respectively, in both the global calibration and the fitting-forecasting study, thereby matching the initial histopathologic assessment available for these tumors. Similar results were obtained for the patients with high-risk prostate cancer in Fig. 9B and Supplementary Fig. S3A. For the latter, our modeling framework also anticipated the progression from GS 4 + 3 to GS 4 + 4 using the personalized predictions at the second imaging timepoint, which was 410 days (i.e., ∼1.1 years) earlier than the final histopathologic assessment at surgery. Importantly, although the early biopsies of the tumor in Fig. 9C report a low-risk case, the personalized biomechanistic model predictions (calculated at the second imaging timepoint) consistently classified this tumor as high risk and the global calibration (calculated at the third imaging timepoint) further confirmed prostate cancer progression right after the second scan. For this patient, our model anticipated the detection of prostate cancer progression at the time of the second mpMRI, which was 1,011 days (i.e., ∼2.8 years) before the final histopathologic assessment after surgery. Similar results are obtained for the prostate cancer case in Supplementary Fig. S3B, for which tumor progression was detected 677 days (i.e., ∼1.9 years) before surgery. Finally, the low-risk patient in Supplementary Fig. S3C was consistently identified as high-risk with the predictions of the fitting-forecasting study, but assimilation of the third mpMRI dataset in the global calibration scenario correctly classified the tumor as low-risk.

The majority of patients with newly diagnosed prostate cancer are estimated to exhibit low or intermediate-risk prostate cancer at diagnosis (5, 6), for which AS is a standard clinical option (2, 3, 12). To optimize patient monitoring during AS and rationalize the treatment of prostate cancer, it is of utmost importance to accurately identify the patients with indolent disease, who can continue to reap the benefits of AS and avoid treatment, as well as the patients with progressing prostate cancer, who may require immediate treatment. Nevertheless, current AS protocols largely rely on assessing prostate cancer status at standardized frequencies that are set upon the average dynamics of the disease observed in large clinical studies (2, 12, 13). This population-based, observational approach cannot anticipate prostate cancer progression and does not enable an optimal surveillance of the unique growth dynamics of each patient's tumor. To address these limitations of AS, we propose to use personalized computational forecasts of prostate cancer growth obtained with a spatiotemporal biomechanistic model informed by the longitudinal mpMRI data collected during AS for each individual patient. Hence, our forecasting technology ultimately aims at advancing AS from its current observational, population-based standard toward a predictive, patient-specific paradigm.

Our prostate cancer forecasts primarily consist of the prediction of spatiotemporal tumor growth in terms of normalized tumor cell density maps over the 3D anatomy of the patient's prostate. In the pilot study presented herein, we have shown that our approach can represent and predict prostate cancer growth according to global metrics of tumor volume and cellularity (see Figs. 5 and 7). This performance is analogous to that reported in prior tumor forecasting studies (18, 28, 29). Our results also show that our computational framework exhibits a promising potential to capture the evolving 3D morphology of the tumor within the patient's prostate (see Figs. 4 and 6). This information can be central to the precise planning of radical treatments, such as surgery or radiotherapy (17, 18, 22, 23, 29). In addition, computational predictions of the 3D tumor morphology within the patient's prostate can be leveraged to guide biopsies and, therefore, reduce the risk of underestimation of GS with respect to the histopathology analysis of the surgically-removed prostate (2, 11). Furthermore, despite the simplicity of the biomechanistic model leveraged in this study, our results show a promising trend in the local agreement between imaging measurements and model calculations of the normalized tumor cell density maps (see Figs. 4–7), which is comparable to other tumor forecasting studies (18, 23, 28, 29, 31).

Our patient-specific forecasts further enable the prediction of prostate cancer risk by leveraging a set of biomarkers calculated from normalized tumor cell density maps obtained from the personalized model. In particular, in this pilot study we found that the mean proliferation activity of the tumor can be used as a promising biomarker of higher-risk prostate cancer, whose classifying performance improved in combination with the total tumor index (see Eqs. (9) and (10)). Tumor cell proliferation has been found to be a central mechanism driving tumor growth and treatment response in previous biomechanistic modeling studies, and model-based biomarkers of proliferation activity have been found to correlate with tumor progression and unfavorable pathologic outcomes (22, 25, 26, 30, 54, 56). In the case of prostate cancer, previous clinical studies have also shown that a high proliferation activity has been associated with higher GS as well as with poorer prognosis, treatment outcomes, and survival (57, 58). Furthermore, the total tumor index combines three key variables in the analysis of prostate cancer growth in the clinical literature: cellularity (measured via the average normalized tumor cell density), tumor volume, and prostate volume (see Eq. (9)). While the first two have been directly associated with higher-risk prostate cancer (4, 12, 13, 41–48), the latter has been suggested to have an inverse correlation (i.e., such that larger prostates tend to harbor tumors with more pathologically favorable features; see refs. 32, 55). Thus, the model-based quantities selected in the current study to construct a predictive classifier of high-risk prostate cancer agree with previous analyses of the driving mechanisms of tumor growth in modeling studies, as well as observations of prostate cancer biology and clinical progression in the literature. These results along with the promising performance of the bivariate classifier to predict progression at the time of the second mpMRI scan suggest that our forecasts of prostate cancer risk could contribute to guide fundamental clinical decisions in AS, such as the frequency of monitoring tests and the best timing to direct a patient with progressing prostate cancer to definitive treatment. For example, should the prostate cancer risk prediction reveal that a tumor will be at low risk at the 1- to 3-year horizon, then the interscan time can be extended to monitor the tumor. Conversely, if the forecast of prostate cancer risk reveals a trend toward progression, the subsequent scans could be performed earlier to confirm this pathologic event and promptly treat the tumor.

Despite the promising results in personalized forecasting of prostate cancer growth during AS observed in this pilot study, we also acknowledge several limitations. First, we employed a small patient cohort (n = 7). Beyond statistical limitations (e.g., potential data correlation and low statistical power), this small cohort is not representative of the vast intratumoral and intertumoral heterogeneity of the newly detected prostate cancer cases that are eligible for AS. In addition, the limited size of the cohort also prevented us from training and validating the prostate cancer risk classifiers in different patient subgroups, constructing classifiers with more than two input model-based biomarkers, and defining a multitiered classifier separately accounting for tumors with GS 3 + 3, 3 + 4, and 4 + 3 or higher. Further work addressing this latter limitation would be of great clinical interest, because some GS 3 + 4 tumors are considered eligible for AS (2, 12). Thus, future prostate cancer forecasting studies should involve a larger cohort with more diversity of tumors, for example, in terms of location, number and size of lesions, prostate volumes, and GS changes during AS (e.g., stable indolent disease, early versus late progression, overall GS values; see refs. 2, 11–15). Importantly, these studies will rely on the computational pipeline presented in this pilot study, which demonstrated that mpMRI-informed model provides reasonable results in predicting prostate cancer growth and progression. Indeed, this approach of performing a pilot study followed by larger population studies has been a successful strategy to develop computational tumor forecasting technologies for breast and brain cancers (18, 20, 21, 28, 29). Of note, our predictions of spatiotemporal prostate cancer growth are obtained from a biomechanistic model that is calibrated on a patient-specific manner. Hence, only the longitudinal data from each individual patient is needed to constrain the model and make personalized forecasts, whose validation only depends on the subsequent mpMRI data collected from the patient during the course of AS. While a large cohort may provide more information of our biomechanistic model's predictive performance over more diverse prostate cancer cases, the calculation and validation of additional forecasts of spatiotemporal growth of prostate cancer would still be patient-specific (i.e., our biomechanistic model does not necessitate a cohort of patients for training, contrary to pure data-driven approaches (27, 59–61)). Conversely, our predictions of prostate cancer progression rely on a data-driven approach (i.e., logistic classifiers) and, given the small size of the patient cohort used in our study, we did not define separate training and validation cohorts. Despite this limitation, the promising preliminary results on prostate cancer progression presented in this work now justify performing subsequent studies in large patient populations featuring different training and validation cohorts to assess the performance, reproducibility, and reliability of our personalized forecasts and model-based biomarkers. Second, we used both biopsy and surgically-based measurements of GS to train our high-risk prostate cancer classifiers. However, prostate cancer biopsies may underestimate GS mainly due to sampling bias and the incomplete understanding of imaging-histopathology data correlation (3, 11–13, 16, 48). Thus, a more robust training of our classifiers could be achieved by exclusively leveraging GS measurements on surgically removed prostate specimens paired to recent presurgery mpMRI scans.

Another limitation is that our personalized forecasting approach does not account for uncertainty quantification and requires at least two mpMRI datasets to initialize and calibrate the model. For the pilot study presented in this work, we employed a deterministic approach following other previous studies pioneering the use of tumor forecasting in clinical settings (18, 21, 22, 28–31). Hence, it suffices to work with three mpMRI datasets per patient. This amount of data enables the analysis of our tumor forecasting method in a global calibration scenario (in which we assess whether we can reproduce observed tumor growth and progression using all available mpMRI data) as well as in a fitting-forecasting scenario (in which we assess whether we can predict tumor growth and progression at the time of the last mpMRI scan by informing the model with only the first two mpMRI datasets). We believe that the results obtained in these two scenarios provide a good initial assessment of our patient-specific tumor forecasting methods and motivate further development and validation over larger cohorts. Of note, despite the diverse interscan times during AS for each patient (which ranged from 0.5 to 2.6 years), we would like to remark that our methods show strong predictive capabilities regardless of the time between the mpMRI datasets used for calibration as well as of the time between the last mpMRI dataset used for the personalized calibration of our biomechanistic model horizon and the forecasting horizon (i.e., the times of the second and third mpMRI datasets for each patient in this study, respectively). Nevertheless, to maximize the utility of our predictive technology during AS, it would be necessary to enable tumor forecasting from the first imaging timepoint and then update the model as more data become available. Indeed, the comparison of the results of the global calibration study (three mpMRI scans per patient) and those from the fitting-forecasting study (two imaging datasets per patient) suggests that incoming longitudinal data enable the refinement of the model parameterization and a more accurate representation of tumor morphology closer to the calibration horizon (e.g., lower RMSE or larger DSC, as observed herein). We believe that a probabilistic framework (18, 23) would enable a robust implementation of our predictive technology including the uncertainty quantification of input data and model outcomes (e.g., parameters, model-based biomarkers), the update of the distributions of model parameters and predictions with new data collected during the course of AS, and tumor forecasting from the first imaging session. In such a probabilistic framework, the model provides a collection of tumor forecasts for multiple samples of the distribution of the model parameters based on existing data for each patient and prior knowledge over the population, thereby enabling predictions from baseline. Then, these forecasts can be analyzed using standard statistical metrics (e.g., median and range, mean and 95% confidence interval) and risk-based metrics (e.g., quantiles and superquantiles) that enable uncertainty quantification. In particular, the latter can inform about the risk of tumor cell density or other tumor biomarker reaching thresholds suggesting tumor progression (62). Furthermore, as more longitudinal data are collected for the patient, the parameter distributions can be updated, and this additional information can be transferred to update the distributions of tumor forecasts and biomarkers via Bayesian data assimilation (27, 62, 63). Because reproducibility and reliability are substantial barriers for the clinical adoption of new imaging biomarkers, uncertainty quantification could also contribute to the clinical translation of the model-based biomarkers proposed herein under variations in image acquisition as well as segmentation and registration errors.

Our model is only informed by mpMRI data, which constitutes another limitation of this work. Despite the wealth of spatiotemporal information provided by this imaging modality, patients with prostate cancer in AS are also monitored using PSA tests, which are often performed more frequently than mpMRI scans (2, 12, 13). Thus, combining the spatiotemporal model of tumor cell dynamics used here (i.e., Eqs. (2) or (3)) with an equation describing PSA dynamics could facilitate a more frequent update of the patient-specific parameters and include PSA dynamics in the construction of the prostate cancer risk classifier (17, 32). The reason for not including serum PSA in the current model is 2-fold. First, we wanted to build a parsimonious model that could reliably reproduce and predict spatiotemporal prostate cancer growth, for which mpMRI data suffice, and progression, for which GS data are necessary. Second, serum PSA is a scalar value and tumor cell density is a spatial quantity. Although we have already proposed a method to relate these two concepts (17), further studies are needed to investigate the joint calibration of the parameters of the equations governing the temporal dynamics of serum PSA and the spatiotemporal dynamics of tumor cell density, along with the analysis of their contributions to predict prostate cancer risk. In addition, genomic classifiers are gaining increased interest for risk stratification of newly diagnosed prostate cancer. These tests have been strongly and independently associated with adverse pathology, treatment failure, metastasis, cancer-specific and overall survival, and conversion from AS to definitive treatment (64–66). While their use is not standardized in AS protocols and they are seldom performed longitudinally, the results of these genomic classifiers could be combined with our model-based biomarkers to refine the prediction of prostate cancer progression during AS as well as to refine the construction of prior distributions of parameter values within a probabilistic forecasting framework.

Finally, the simplicity of the model is a central limitation of this study. In particular, the biomechanistic model used herein tends to homogenize the intratumoral spatial distribution of the normalized tumor cell density map. This effect is caused by the Fisher–Kolmogorov equation, which we used to describe untreated prostate cancer growth in our biomechanistic model (see Eqs. (2) and (3)). The natural solution of this equation is a travelling wavefront (51) that volumetrically expands a spatial region with nonzero tumor cell density to occupy other neighboring healthy regions, such that normalized tumor cell density progressively increases in the central part of the tumor and decreases along the travelling wavefront. Hence, although the initial conditions defined upon the normalized tumor cell density map extracted from the first mpMRI of each patient introduce information about the internal heterogeneity of the tumor architecture, the ensuing model simulation will progressively adapt the internal spatial distribution of normalized tumor cell density to match the natural solution of the Fisher–Kolmogorov equation. This effect contributes to explain the comparatively poorer performance of the personalized model at the third imaging timepoint with respect to the second one in both the global calibration and fitting-forecasting studies (see Figs. 4 and 6), as well as the relatively lower local model-data agreement with respect to previous studies employing more advanced models (18, 23, 28, 29, 31). In addition, the progressive model-driven spatial homogenization of the normalized tumor cell density maps precludes accurate predictions of local GS values, for example, using the hyperbolic tangent mapping built to support the imaging preprocessing pipeline (see Fig. 2). The analysis of these local GS maps would also require data on the spatial position of the prostate cancer–positive biopsy cores and surgical specimen regions, but this information was not available for this study. Nevertheless, the prediction of local GS maps would be especially important to precisely identify the regions of the tumor for biopsy (2, 3, 16), to adapt the treating fields during radiotherapy (22, 23, 29), and, in general, to provide richer predictions of prostate cancer progression during AS. Thus, future studies should investigate model extensions that incorporate mpMRI-informed spatial mechanisms that aim at retaining and evolving the intratumoral heterogeneity observed in imaging measurements of normalized tumor cell density maps. To achieve this goal, previous studies have leveraged heterogeneous and anisotropic parameterizations (18, 21, 23, 28, 29, 31), mechanical constraints to tumor cell mobility and proliferation (18, 21, 28, 29, 31), and information on tumor-supporting vascularity (18, 20, 21, 28, 67). In particular, given that the prostate is a confined organ in the pelvic region and that both a tumor and potentially coexisting BPH may enlarge the prostate volume, the mechanical stress field within the organ can have a major impact on tumor dynamics (32). Thus, the extension of the model to include the mechanical effects on prostate cancer growth constitutes a promising avenue of research. Toward this end, it would also be necessary to obtain a segmentation of the central gland of the prostate in T2W images to define the mechanical stress due to BPH (32), which can be straightforwardly processed along the corresponding prostate segmentation within our computational pipeline. Another promising extension to refine the spatial heterogeneity of our prostate cancer forecasts would be the definition of a spatially-varying carrying capacity map (θ(x)) based on local measurements of vascularization obtained via dynamic contrast-enhanced MRI (67) and the local level of mechanical stress, which may limit tumor cell density via changes in prostatic tissue architecture caused by prostate cancer and BPH (4, 18, 32, 46–48, 55, 68, 69).

In the future, we believe that further developing our computational modeling framework addressing the aforementioned limitations will enable the construction of a digital twin (27, 62) for each individual patient enrolling in AS. These technologies not only enable a seamless integration of patient data and personalized model forecasts, but also exploit them to define specific metrics that aid the treating physician in making decisions about the timing and type of further testing or treatments, thereby providing a robust and systematic framework for the optimal management of each patient's tumor. By further extending our biomechanistic model to describe the response to standard first-line treatments for prostate cancer after AS, the digital twin could also be trained to calculate therapeutic outcomes or expected survival, and, eventually, optimal regimens maximizing them while minimizing side effects (18, 20–33, 62). Future technological developments in mpMRI can also contribute to new and more accurate imaging measurements of prostate cancer (3, 70), which can contribute to better inform our model and enhance the reliability of our model predictions. Thus, we posit that the work presented in this study constitutes an important first step toward the construction of a digital twin enabling a personalized, predictive paradigm for the clinical management of prostate cancer.

We have presented a clinical-computational framework for the personalized forecast of untreated, newly diagnosed prostate cancer growth during AS constructed upon a spatiotemporal biomechanistic model informed by patient-specific longitudinal mpMRI data. This technology tracks and predicts the 3D development of the tumor within the patient's prostate geometry, while also enabling to anticipate the progression to higher risk-disease. Although further model development and investigation within larger cohorts are necessary, we posit that the personalized predictive technology presented in this pilot study can contribute to guide optimal clinical decision-making regarding the frequency and type of monitoring during AS, as well as the optimal time, modality, and planning of ensuing treatment.

G. Lorenzo reports grants from European Commission during the conduct of the study; in addition, G. Lorenzo has a patent to US20230274842A1 pending, a patent to US20190198177A1 issued, and a patent to WO/2018/037137 pending. J.S. Heiselman reports grants from NIH during the conduct of the study. H. Gomez reports grants from National Science Foundation during the conduct of the study; grants from Eli Lilly outside the submitted work; in addition, H. Gomez has a patent to US Patent Application Number 16327875, with publication number US20190198177A1 issued, a patent to US Patent Application Number 18135580, with publication number US20230274842A1 pending, and a patent to International patent WO/2018/037137 pending. T.E. Yankeelov reports grants from Cancer Prevention and Research Institute of Texas during the conduct of the study; in addition, T.E. Yankeelov has a patent to US20230274842A1 pending. A. Reali reports a patent to US20230274842A1 pending. T.J.R. Hughes reports a patent to WO/2018/037137 pending, a patent to US20190198177A1 issued, and a patent to US20230274842A1 pending. No disclosures were reported by the other authors.

G. Lorenzo: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing-original draft, project administration, writing-review and editing. J.S. Heiselman: Resources, data curation, software, formal analysis, validation, investigation, methodology, writing-original draft, writing-review and editing. M.A. Liss: Resources, data curation, investigation, methodology, writing-original draft, writing-review and editing. M.I. Miga: Formal analysis, investigation, methodology, writing-review and editing. H. Gomez: Conceptualization, formal analysis, investigation, methodology, writing-review and editing. T.E. Yankeelov: Conceptualization, formal analysis, investigation, methodology, writing-review and editing. A. Reali: Conceptualization, formal analysis, funding acquisition, investigation, methodology, writing-review and editing. T.J.R. Hughes: Conceptualization, formal analysis, investigation, methodology, writing-review and editing.

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 838786. J.S. Heiselman and M.I. Miga thank the NIH for funding through R01EB027498. H. Gomez was partially funded by the National Science Foundation (CMMI 1852285). We also thank the Cancer Prevention and Research Institute of Texas (CPRIT) for support through CPRIT RR160005; T.E. Yankeelov is a CPRIT Scholar in Cancer Research. A. Reali acknowledges the support of the Italian Ministry of University and Research (MUR) through the PRIN project COSMIC (No. 2022A79M75), funded by the European Union – Next Generation EU. We thank the Texas Advanced Computing Center (TACC) for providing high-performance computational resources that contributed to the results presented in this work. The opinions, findings, and conclusions, or recommendations expressed are those of the authors and do not necessarily reflect the views of the funding agencies.