Modeled Early Longitudinal PSA Kinetics Prognostic Value in Rising PSA Prostate Cancer Patients after Local Therapy Treated with ADT +/− Docetaxel

Simple Summary In prostate cancer patients with rising prostate-specific antigens (PSAs) after primary local therapy, and at high risk of metastatic disease based on a high Gleason score or ISUP grades 4–5 and/or short PSA doubling time, androgen-deprivation therapy (ADT) was shown to be effective with or without salvage radiotherapy. In the present paper, the analysis of a phase III trial dataset comparing ADT +/− docetaxel demonstrates that longitudinal PSA kinetics during the first 100 days of treatment, assessed using mathematical modeling, is associated with patient prognosis. Indeed, multivariate analyses showed that the modeled PSA production rate constant KPROD, and the elimination rate constant KELIM, exhibited strong and independent prognostic values regarding PSA progression-free survival (PSA-PFS) and overall survival (OS), respectively. As it was shown with the longitudinal CA-125 kinetic in ovarian cancers, the modeled longitudinal PSA kinetics during the 100 days of treatment may represent a novel prognostic factor in patients with metastatic prostate cancers. Abstract Background: In metastatic prostate cancer (PCa) patients, androgen-deprivation therapy (ADT) combined with chemotherapy or next-generation androgen receptor targeted agents is a new standard treatment. The objective of the present study is to assess longitudinal PSA kinetics during treatment using mathematical modeling, to identify the modeled PSA kinetic parameters able to exhibit early prognostic/predictive values. Methods: Phase III clinical trial dataset (NCT00764166) comparing ADT +/− docetaxel in 250 locally treated patients for PCa with rising PSA levels, who were at high risk of metastatic disease was assessed. A kinetic-pharmacodynamic (K-PD) model was used to fit PSA kinetics during the first 100 treatment days, to estimate the modeled PSA production rate K (KPROD) and elimination constant rate K (KELIM). The prognostic value of these parameters, considered as categorized (favorable vs. unfavorable) covariates regarding PSA progression-free survival (PSA-PFS) and overall survival (OS), was assessed using univariate/multivariate analyses. Results: Data from 177/250 patients was assessed. KELIM exhibited a significant prognostic value regarding PSA-PFS and KPROD regarding OS (univariate analysis). In the PSA-PFS final multivariate model, KELIM and the primary therapy type were significant. The OS multivariate model integrated both KPROD and baseline PSA doubling-time. Conclusion: In this first study assessing the modeled PSA kinetics prognostic value in PCa patients treated with systemic treatments, KELIM and KPROD exhibited respective prognostic values regarding PSA-PFS and OS.


Introduction
In patients treated with primary therapy (surgery or radiotherapy) for localized prostate cancer (PCa), biochemical relapses (BR) with the rise in the isolated prostatespecific antigen (PSA) are observed in up to 30% of them within 7 years following the primary therapy [1]. For men who underwent radical prostatectomy (RP) as the primary treatment and subsequently developed PSA recurrence, the main unfavorable prognostic factors were a PSA doubling time (PSA-DT) of 1 year and a pathological Gleason score (pGS) of 8-10 (International Society of Urological Pathology (ISUP) grades [4][5]. For patients with a PSA recurrence following primary RP, an interval from primary therapy to biochemical failure (IBF) of <18 months and a biopsy Gleason score (bGS) of 8-10 (ISUP grades [4][5] were the main unfavorable prognostic factors. The EAU prostate cancer guidelines panel recommended the use of a novel BR classification system that stratified patients with BR into low-vs. high-risk categories. After BR in case of RP, salvage radiation therapy can be an option in high-risk patients [2]. For decades, androgen-deprivation therapy (ADT) was considered as a standard treatment in the case of metastatic hormone-sensitive prostate cancer (mHSPC). Over the last 5 years, the treatment landscape has dramatically changed with the addition of 4 systemic agents that previously demonstrated benefits in the castrate-resistant setting [3][4][5][6][7]. In patients with only a biochemical relapse, no systemic treatment was approved [8]. However, the high risk of developing metastases may justify early ADT administration [9], which can be balanced with the long term tolerability of this treatment [10].
In metastatic hormone-naive prostate cancer patients, docetaxel and abiraterone became the treatment options, based on CHAARTED (Chemo Hormonal Therapy versus Androgen Ablation Randomized Trial in Extensive Disease), LATITUDE and STAMPEDE (Systemic Therapy in Advanced and Metastatic Prostate Cancer evaluation of Drug Efficacy) [11][12][13][14]. In addition, recently enzalutamide and apalutamide were added to the armamentarium, after ARCHES, ENZAMET and TITAN phase III studies [15][16][17]. Radiation treatment of the primary in patients with oligometastatic disease recently became an additional option for patients, due to STAMPEDE-arm H. In line with the data suggesting the overall survival (OS) benefit with the addition of docetaxel to ADT in patients with metastatic hormone-naive prostate cancer (GETUG-15, STAMPEDE and CHAARTED trials) [11,12,18], Oudard et al. investigated the utility of this strategy in patients with rising PSA, who are hormone-naive at a high-risk of metastatic progression in a randomized trial. Unfortunately, no benefit in PSA-progression free survival (PSA-PFS) and OS was found with the addition of docetaxel to ADT [19].
Despite this outcome, it is important to search for subgroups of patients who may have benefited from this approach. Consistent with the strategy developed with the modeled kinetics of the serum tumor marker CA-125 in ovarian cancer patients treated with chemotherapy, producing a strong and reliable prognostic marker of OS called KELIM (CA-125 ELIMination constant rate K) [20][21][22], we proposed to assess the prognostic value of PSA kinetics in this trial as a way of identifying some patients who might have benefited from adding docetaxel to ADT.
Of note, the PSA working group has adopted a PSA response criterion, defined as a >50% decrease in PSA with an absolute PSA level lower than 0.2 ng/mL [23]. However, as already found with CA-125 in ovarian cancers, this two time point-based kinetic strategy is largely impacted by the high inter-and intra-individual variability of sparse PSA concentrations measured with different assays, far from the actual dynamics of serum marker kinetics, thereby limiting the comparability in different patient populations than those assessed initially [24].
In that context, mathematical modeling with population kinetic approach of serum tumor marker dynamics during treatment offers several advantages, as it enables to calculate the mathematical equations describing the longitudinal serum tumor kinetics of the population of patients, and then those of individual subjects based on a minimum of three values, with reduced impact of measurement variability. Moreover, it enables to extract a modeled kinetic parameter of interest prone to exhibit prognostic or predictive value, as seen with CA-125 KELIM [20][21][22].
The purpose of the present study is to assess the longitudinal PSA kinetics using mathematical modeling in patients enrolled in the Oudard et al. trial, in order to identify the modeled kinetic parameter of interest (such as KELIM) prone to exhibit prognostic value or predictive value regarding the benefit from adding docetaxel to ADT, in terms of PSA-PFS and OS.

Patients
We analyzed the dataset of the randomized phase III trial (NCT00764166) [19], comparing two arms of treatment (endocrine therapy with ADT, or ADT + docetaxel 70 mg/m 2 every 3 weeks) in 250 patients with localized PCa previously treated with primary radical prostatectomy (RP) or radiotherapy (RT) and experiencing isolated PSA rising (at least 3 consecutives PSA measurements, such as defined by the American Society for Therapeutic Radiology and Oncology PSA recurrence criteria [25]) with at least one or more high-risk of metastatic disease characteristics (Gleason score ≥ 8, PSA velocity > 0.75 ng/mL, PSA doubling-time (DT) < 6 months, node-positive adenocarcinoma, positive surgical margins, and time to PSA recurrence 12 months or less) [26]. The primary endpoint of the trial was PSA-PFS, while secondary endpoints were PSA response, radiologic PFS, and OS.
To be included in this analysis, a minimum of three PSA assays in the first 100 days of treatment was required to be able to fit individual PSA kinetics to a non-linear model. The patients who had PSA titers rising >10% during the first 100 days of treatment were excluded, as this treatment could not be considered effective, which was inadequate with respect to the objective of the present study.
In addition to patients PSA concentrations and sampling dates, the following clinical data were collected: age at diagnostic, treatment arm, baseline PSA doubling-time (PSA-DT), baseline PSA velocity, the type of primary therapy (RT/RP), pTNM stage, and Gleason score [19].

Mathematical Modeling of PSA Kinetics
Individual PSA data were log transformed to eliminate right-skewness in the data distribution. The mathematical modeling of PSA kinetics in the blood with a non-linear mixed-effect model was the same used from CA-125 in ovarian cancer [21,22], and assessed with the observed PSA values measured during the first 100 days of treatment starting on cycle 1 day 1 of chemotherapy, regardless of treatment arm. PSA kinetics were described by the PSA production rate KPROD (ng·mL −1 ·day −1) , a PSA elimination constant rate KELIM (day −1 ), both being impacted by the indirect effects of treatment with ADT +/− docetaxel ( Figure 1).
The semi-mechanistic model structure relied on a central compartment where the blood serum concentration of PSA was described through a production rate KPROD, balanced with an elimination rate constant KELIM. KPROD was regulated by the effects of the systemic treatment, described with a 2-virtual-compartment model (a central compartment C1 receiving chemotherapy dosing and a transit compartment C2) on the cancer cells through an indirect effect, characterized by the half-maximal effective concentration EC50. Because the concentrations of the chemotherapy were not available, a kineticpharmacodynamic (K-PD) model was developed to assess the longitudinal kinetics of PSA with the administration dose of chemotherapy arbitrarily set to 1, as already performed for the pharmacokinetics model [27]. The semi-mechanistic model structure relied on a central compartment where the blood serum concentration of PSA was described through a production rate KPROD, balanced with an elimination rate constant KELIM. KPROD was regulated by the effects of the systemic treatment, described with a 2-virtual-compartment model (a central compartment C1 receiving chemotherapy dosing and a transit compartment C2) on the cancer cells through an indirect effect, characterized by the half-maximal effective concentration EC50. Because the concentrations of the chemotherapy were not available, a kinetic-pharmacodynamic (K-PD) model was developed to assess the longitudinal kinetics of PSA with the administration dose of chemotherapy arbitrarily set to 1, as already performed for the pharmacokinetics model [27].

Model Qualification
The PSA predictions, obtained with the fitting of data on the K-PD model, were qualified using common algorithms, such as relative standard error (RSE), goodness-of-fit plots (GOF), meaning plots comparing observations vs. predictions, the distribution of normalized prediction distribution errors (NPDE) or visual predictive check (VPC). VPC represented 500 replicates of all individual PSA decline profiles simulated from the final parameter estimates of the model: the 10th, 50th, and 90th percentiles of the observed PSA values were visually included in the 95% Confidence Interval (CI) of simulations. VPC were used as models for internal validation.

Prognostic Value of Modeled PSA Kinetics During the First 100 Days
The prognostic value of modeled PSA kinetic parameters regarding PSA-PFS and OS was assessed using univariate and multivariate survival tests. The modeled PSA kinetic parameters of interest (KPROD and KELIM) were standardized and dichotomized by their median, respectively.
Kaplan-Meier and log-rank tests were used for univariate analyses. All variables found significant with p < 0.2 in univariate analyses were integrated in a multivariate Cox proportional hazards regression model. The final Cox regression survival model was obtained using a backward selection. The following potential prognostic factors were also tested in the multivariate analysis: age at diagnostic, treatment arm, baseline PSA doubling time (PSA-DT), baseline PSA velocity, type of previous primary therapy (RT/RP), pTNM stage, and ISUP classification.

Model Qualification
The PSA predictions, obtained with the fitting of data on the K-PD model, were qualified using common algorithms, such as relative standard error (RSE), goodness-of-fit plots (GOF), meaning plots comparing observations vs. predictions, the distribution of normalized prediction distribution errors (NPDE) or visual predictive check (VPC). VPC represented 500 replicates of all individual PSA decline profiles simulated from the final parameter estimates of the model: the 10th, 50th, and 90th percentiles of the observed PSA values were visually included in the 95% Confidence Interval (CI) of simulations. VPC were used as models for internal validation.

Prognostic Value of Modeled PSA Kinetics during the First 100 Days
The prognostic value of modeled PSA kinetic parameters regarding PSA-PFS and OS was assessed using univariate and multivariate survival tests. The modeled PSA kinetic parameters of interest (KPROD and KELIM) were standardized and dichotomized by their median, respectively.
Kaplan-Meier and log-rank tests were used for univariate analyses. All variables found significant with p < 0.2 in univariate analyses were integrated in a multivariate Cox proportional hazards regression model. The final Cox regression survival model was obtained using a backward selection. The following potential prognostic factors were also tested in the multivariate analysis: age at diagnostic, treatment arm, baseline PSA doubling time (PSA-DT), baseline PSA velocity, type of previous primary therapy (RT/RP), pTNM stage, and ISUP classification.
Moreover, concordance index (C-index) as a function of the time was used to assess the prediction benefit related to their addition in the multivariate model.
All survival analyses were performed with a 100-day landmark time point analysis, meaning that all patients experiencing events (lost to follow-up, progression or death for PSA-PFS; and lost to follow-up or death for OS) during the first 100 days were excluded. Indeed, the modeled kinetic parameters were estimated during the first 100 days of treatment. This strategy avoided any bias between survival predictions and PSA kinetics parameters estimations [21,22].

Statistical Analyses and Computing Process
All tests were implemented using a two-sided 0.05 alpha risk. The NONMEM 7.5.0 (ICON Development Solutions, Ellicott City, MD, USA) software was used to fit the semimechanistic model to PSA kinetics data [28]. Parameters were estimated adjusting data to the model with the First Order Conditional Estimates with Interaction (FOCEI) algorithm plus the Laplacian option. Error model included the PSA Lower Limit of Quantification (LLOQ = 0.1 ng/mL) thanks to the M3 method [29]. Standard Errors (SE) were estimated by 2000 samples bootstrap. Survival analyses and graphical representations of results were performed in R language (R© 3.6.1 software). The Xpose 4.7.1 package in R was used for representation of model qualification (visual predictive check (VPC) and goodness-of-fit (GOF) plots) [30].

Patient Selection
Of the 250 patients enrolled in the trial, 177 patients were assessable for PSA kinetic modeling, including 94 and 83 patients in ADT + docetaxel and ADT alone arms, respectively. The characteristics of the assessed patients are presented in Supplementary Table S1. In summary, a bit more than 50% of patients had previous T3/T4 stage diseases. Most of them were previously treated with prior prostatectomy (68%), and/or radiotherapy (31%). Moreover, the majority of patients had high-risk features at inclusion, based on PSA doubling time and PSA velocity (Supplementary Table S1).
The characteristics of patients were well balanced in between treatments arms. The features of excluded patients were similar to those assessed in the present study for most of them, except for baseline PSA and percentage of patients previously treated with radiotherapy, which were slightly lower in the excluded patients.
The median follow-up was 9.8 years, and the median number of PSA values by patients was 4. Of the 177 patients eligible for the PSA kinetics modeling, 173 were assessable for PFS survival analyses, according to the 100-day landmark, including 92 and 81 patients in ADT + docetaxel and ADT alone arms, respectively. No patient presented an OS event before 100 days (Figure 2).

Model Qualification
Parameter estimates from the semi-mechanistic model are presented in the Supplementary Table S2. Relative standard errors of parameters and their inter-individual variation, which represent the estimation precision, were lower than 25%, as was the shrinkage (<30%, excepted for K and EC50), thereby suggesting limited risks of biased individual estimates of parameters.
The GOF plots showed that individual PSA profiles were properly fitting observations during the first 100 treatment days by the model. Weighted residuals approximately follow a normal distribution. Moreover, the VPC revealed the correct predictive performance of the model with no aberrant predictions (Supplementary Figure S1).

KELIM and KPROD Prognostic Values on PSA-PFS and OS in Univariate Analyses
Among the five parameters estimated, KELIM and KPROD prognostic values were considered of interest, as these parameters, meant to characterize PSA production and elimination, were well estimated compared to other modeled kinetic parameters, with lower RSE and lower shrinkage (Supplementary Table S2). KPROD and KELIM were standardized by their median, respectively (median KPROD = 0.002 ng·mL −1 ·day −1 , median KELIM = 0.081 day −1 ), and only the standardized parameters were used in the following analyses.

Model Qualification
Parameter estimates from the semi-mechanistic model are presented in the Supplementary Table S2. Relative standard errors of parameters and their inter-individual variation, which represent the estimation precision, were lower than 25%, as was the shrinkage (<30%, excepted for K and EC50), thereby suggesting limited risks of biased individual estimates of parameters.
The GOF plots showed that individual PSA profiles were properly fitting observations during the first 100 treatment days by the model. Weighted residuals approximately follow a normal distribution. Moreover, the VPC revealed the correct predictive performance of the model with no aberrant predictions (Supplementary Figure S1).

KELIM and KPROD Prognostic Values on PSA-PFS and OS in Univariate Analyses
Among the five parameters estimated, KELIM and KPROD prognostic values were considered of interest, as these parameters, meant to characterize PSA production and elimination, were well estimated compared to other modeled kinetic parameters, with lower RSE and lower shrinkage (Supplementary Table S2). KPROD and KELIM were standardized by their median, respectively (median KPROD = 0.002 ng·ml -1 ·day -1 , median KELIM = 0.081 day −1 ), and only the standardized parameters were used in the following analyses.
KELIM and KPROD estimations were not significantly different across treatment arms (KELIM, p = 0.39; KPROD, p = 0.093) (Figure 3A,B).   (Figure 4A,B). KELIM was also significant for OS, with a median of 157 months in patients with unfavorable KELIM, compared to a median survival "not reached" in patients with favorable KELIM (p = 0.035) ( Table 1). KELIM exhibited a significant prognostic value regarding patient PSA-PFS. The median PSA-PFS was 10.1 months in patients with unfavorable KELIM vs. 15.1 months in patients with favorable KELIM (p < 0.001) ( Figure 4A,B). KELIM was also significant for OS, with a median of 157 months in patients with unfavorable KELIM, compared to a median survival "not reached" in patients with favorable KELIM (p = 0.035) ( Table 1).    KPROD also exhibited significant prognostic value regarding patient PSA-PFS and OS ( Figure 4C,D). The median PSA-PFS was 9.8 months for patients with unfavorable KPROD vs. 15.1 months for patients with favorable KPROD (p < 0.001). The median OS was 157 months in unfavorable KPROD patients, and not reached in those with a favorable KPROD (p = 0.037) ( Table 1).

C-Index Analyses
Looking at the C-index, the combination of KELIM and primary therapy seemed to be the most informative model, rather than KELIM alone or primary therapy alone for the PSA-PFS ( Figure 4E).
Regarding OS, the combination of KPROD with baseline PSA-DT was associated with the highest predictive index ( Figure 4F).

Discussion
As for ovarian cancer patients, K-PD modeling was able to characterize the longitudinal PSA dynamics during treatment based on PSA production KPROD, PSA elimination KELIM, and indirect treatment effect on patients, independent of the selected time points used to calculate it. A previous study already suggested that mathematical modeling was a promising tool to assess the longitudinal kinetics of PSA after radical prostatectomy and to predict the risk of early relapse [31].
A major outcome of the present study is the demonstration that modeled PSA kinetics during the first 100 days of systemic treatment inform on the patient prognosis in terms of PSA-PFS and OS. Indeed, patients with favorable KELIM and KPROD experienced significant longer median PSA-PFS and OS in both arms. PSA KELIM was found to be associated with PSA-PFS, but exhibited lower prognostic value regarding OS (contrary to what was reported with CA-125 KELIM in ovarian cancer patients) [21,22]. On the other hand, KPROD seemed to be a stronger prognostic factor regarding OS. Multivariate analyses confirmed that KELIM and KPROD prognostic values regarding PFS and OS were independent, regardless of other prognostic covariates.
Therefore, the present study suggests that KELIM and KPROD might help to identify patients likely to experience shorter or longer survivals.
As KELIM and KPROD values were similar in both treatment arms, these parameters could be calculated in all patients regardless of their current treatment. On the other hand, it could not be used to identify the patients who would benefit from adding docetaxel to ADT. Therefore, it was not a predictive marker of docetaxel efficacy.
Similar outcomes were obtained with the CA-125 modeled kinetic elimination rate constant K (KELIM) in ovarian cancer patients treated with chemotherapy. Indeed, a meta-analysis database of 8 randomized trials involving more than 5000 patients was investigated. KELIM difference between treatment arms was not a surrogate marker of survival benefit in a meta-analysis with more than 5000 patients (Corbaux et al. Proc ESMO 2021). However, consistent with other 7 studies in more than 7000 patients, KELIM was found to be a strong and reproducible prognostic marker regarding PFS and OS, which could be used to discriminate 3 prognostic groups. The poor prognosis group was identified as the subgroup that should be prioritized for innovative treatment meant to reverse the chemoresistance.
As a consequence, additional retrospective analyses on other independent trial datasets will be required to fully understand the prognostic/predictive role of modeled longitudinal PSA kinetic parameters, such as KELIM or KPROD, in patients treated with chemotherapy and/or endocrine therapy, along with prospective validation studies. If the prognostic/predictive value of these parameters was confirmed, they could be easily calculated online at patient bedside on https://www.biomarker-kinetics.org/presentation (accessed on 5 December 2021), where the models can be implemented for routine application.
The limitation of this study is the low number of patients in the trial. This number was reduced further after the selection of patients who experienced a PSA decline in order to obtain a homogeneous population. This strategy was chosen for preventing wrong parameter estimations. However, it might have introduced disbalances in some characteristics of patients and introduced biases in the outcomes.
Cox regression power for OS analysis was 0.38. Our analyses estimated that the number of patients should be doubled in both KELIM/KPROD groups to ensure a standard power of 80%. Thus, our outcomes are mainly hypothesis generating, and the validation of these present results on other PSA databases with larger numbers of patients could reinforce these results.

Conclusions
This is the first study assessing modeled PSA kinetics in recurrent PCa patients with a rising PSA hormone-naive and at high risk of developing metastasis treated with ADT +/− docetaxel.
Despite limitations due to the low number of patients, the outcomes of the present exploratory retrospective study suggest that modeled kinetic parameters KPROD and KELIM might help to identify recurrent PCa patients with a higher probability of biological progression and survival when treated with ADT +/− docetaxel.