Using medical claims database to develop a population disease progression model for leuprorelin-treated subjects with hormone-sensitive prostate cancer.

Androgen deprivation therapy (ADT) is a widely used treatment for patients with hormone-sensitive prostate cancer (PCa). However, duration of treatment response varies, and most patients eventually experience disease progression despite treatment. Leuprorelin is a luteinizing hormone-releasing hormone (LHRH) agonist, a commonly used form of ADT. Prostate-specific antigen (PSA) is a biomarker for monitoring disease progression and predicting treatment response and survival in PCa. However, time-dependent profile of tumor regression and growth in patients with hormone-sensitive PCa on ADT has never been fully characterized. In this analysis, nationwide medical claims database provided by Humana from 2007 to 2011 was used to construct a population-based disease progression model for patients with hormone-sensitive PCa on leuprorelin. Data were analyzed by nonlinear mixed effects modeling utilizing Monte Carlo Parametric Expectation Maximization (MCPEM) method in NONMEM. Covariate selection was performed using a modified Wald's approximation method with backward elimination (WAM-BE) proposed by our group. 1113 PSA observations from 264 subjects with malignant PCa were used for model development. PSA kinetics were well described by the final covariate model. Model parameters were well estimated, but large between-patient variability was observed. Hemoglobin significantly affected proportion of drug-resistant cells in the original tumor, while baseline PSA and antiandrogen use significantly affected treatment effect on drug-sensitive PCa cells (Ds). Population estimate of Ds was 3.78 x 10-2 day-1. Population estimates of growth rates for drug-sensitive (Gs) and drug-resistant PCa cells (GR) were 1.96 x 10-3 and 6.54 x 10-4 day-1, corresponding to a PSA doubling time of 354 and 1060 days, respectively. Proportion of the original PCa cells inherently resistant to treatment was estimated to be 1.94%. Application of population-based disease progression model to clinical data allowed characterization of tumor resistant patterns and growth/regression rates that enhances our understanding of how PCa responds to ADT.

Introduction

Prostate cancer (PCa) is one of the most common cancers among men in the United States, accounting for 44% of all cancer cases along with lung and bronchus, and colorectal cancers [1]. In 2016, there were approximately 180,890 new cases of PCa and 26,120 deaths due to PCa in the United States [1]. Clinically localized PCa is most commonly managed by observation, radical prostatectomy, and radiotherapy (with or without androgen deprivation therapy, or ADT) [2]. However, many patients diagnosed with localized disease ultimately undergo biochemical progression as demonstrated by increasing levels in prostate-specific antigen (PSA), who may then be treated by ADT. ADT acts by depleting gonadal testosterone [3], and it can be achieved by medical castration through the use of luteinizing hormone-releasing hormone (LHRH) agonists or antagonists, or surgical castration (bilateral orchiectomy), both of which are considered equally effective. On the other hand, for patients who present with metastatic PCa, ADT is considered the gold standard of initial therapy that has been shown to reduce tumor-related events [4, 5]. ADT is a widely used systemic treatment for patients with hormone-sensitive PCa, including those who previously receive local therapy as the primary treatment modality and those who present with de novo metastatic disease. However, the majority of patients eventually progress to a castration-resistant state, and it was estimated that the response to ADT typically lasts 14 to 20 months in metastatic PCa [5, 6], and the duration of response to ADT is variable among patients [7]. A few studies have examined clinical factors that predict the time to castration resistance for patients on ADT [8-12], but the detailed time-dependent profile of tumor regression and/or growth in patients with hormone-sensitive PCa on medical ADT has never been fully characterized.

Quantitative population-based disease progression modeling uses mathematical functions and statistical models to describe quantitatively the time course of disease progression in individuals and the entire patient population with or without drug treatment [13, 14]. This modeling approach can trace disease progression over time, and quantify the effects of drug treatment and disease- and patient-specific factors on disease progression to optimize clinical trial design and guide personalized treatment strategies [13]. For example, quantitative disease progression modeling has been applied to clinical data to allow estimation of tumor growth and regression rates in metastatic castration-resistant PCa [14]. Most of the published quantitative population-based disease progression models are developed using data from randomized controlled studies [13-15]. While randomized controlled trials are the golden standard for evaluating new drug treatment and the primary source of clinical research information, they are costly, labor-intensive and time-consuming [16]. In addition, strict inclusion and exclusion criteria are typically used to select subjects in randomized clinical trials within a pre-defined study period, and thus the study results may have limited utility in answering clinical questions in a real-life population [16, 17]. These limitations of randomized controlled trials have led to an increased reliance on using medical claims data in designing healthcare policies related to real-life clinical practice [16, 17]. Medical claims data are a rich and inexpensive source of scientific information to study how a disease responds to medical interventions of interest, and can provide an opportunity to potentially follow a large number of patients for extended periods of time without suffering from high attrition rates [16-18]. To our knowledge, medical claims data have never been used to develop quantitative population-based disease progression models, in which effects of drug treatment and patient-specific factors on disease progression are evaluated and quantified.

The primary goal of this study was to use a medical claims database to construct a quantitative population-based disease progression model for patients with hormone-sensitive PCa on ADT with LHRH agonists. Leuprorelin, the most commonly prescribed LHRH agonist in the Medicare population [19], was selected as the treatment agent of interest in this analysis. Prostate specific antigen (PSA) is an androgen-regulated serine protease produced almost exclusively by cells in the prostate gland [20, 21]. PSA is commonly assessed for PCa screening, and it remains an important marker for monitoring clinical response to ADT treatment [21-24]. Therefore, PSA was used as a molecular marker of tumor burden to measure disease progression of PCa in this study.

Materials and methods

Data source

This study was conducted in commercially insured patients using health claims data provided by Humana covering the period from January 1, 2007 to December 31, 2011. This nationwide database captures anonymized longitudinal, individual-level data on patient demographics, healthcare utilization, inpatient and outpatient diagnostic and procedural codes, laboratory test results, and pharmacy dispensing records for more than 8.1 million commercially insured people in the United States. This study was approved by the Institutional Review Board at the University of Kentucky, and the requirement to document informed consent was waived. All data used in this study were anonymized and deidentified. Data retrieved from SQL queries were pre-processed in R software (version 3.4.3, The R Foundation for Statistical Computing, Austria) for dataset merging, wrangling and formatting (packages used: ‘readr’, ‘dplyr’, ‘stringr’, ‘ggplot2’, ‘data.table’).

Study subjects

Subjects were selected based on the following criteria: 1) diagnosis with malignant PCa (ICD-9-CM code 185 or ICD-10-CM code C61), 2) continuous or intermittent use of leuprorelin as the only agent among LHRH agonists, and 3) availability of a PSA level before initiation of leuprorelin (i.e, a baseline PSA level) and at least one level during treatment.

The following patients were excluded due to inability of our disease progression model to describe certain unusual PSA profiles: 1) PSA levels decreased over time prior to initiation of leuprorelin treatment, the baseline PSA level was undetectable, or there were multiple measurements reported on the same day; 2) subjects had undetectable PSA levels throughout leuprorelin treatment. Patients were also excluded if 3) they had incomplete demographic data (age, region and race). Lastly, patients were excluded if 4) they had extremely low hemoglobin levels (< 6 g/dL), as acute illness was likely involved (S1 Fig).

For all selected subjects, PSA levels starting from baseline to the last available measurement or before the initiation of continuous antiandrogen therapy, surgery, radiotherapy or chemotherapy (whichever occurred first) were included for model development. For each subject, race, age, region and the use of antiandrogens (bicalutamide, enzalutamide, flutamide, nilutamide) within 30 days of leuprorelin initiation (presumably for preventing flare reactions associated with leuprorelin treatment [6]) were extracted for model development. In addition, laboratory measurements including those that have been reported to affect clinical outcomes in castration-resistant PCa were extracted for each patient. These laboratory measurements included aspartate transaminase (AST), alanine transaminase (ALT), serum creatinine (SCR), alkaline phosphatase (ALP), albumin (ALB) and hemoglobin (HGB) [8, 9, 25–33].

Structural model

All modeling steps were conducted in NONMEM (ver. 7.3, ICON Development Solutions, USA). Two mathematical models with different resistance development patterns were developed to describe the observed PSA kinetics. In the first model (Model I), drug resistance was developed from PCa tumor cells that were initially sensitive to LHRH treatment and caused by adaptive responses, such as target gene mutations, altered expression of therapeutic targets, and stimulation of compensatory signaling pathways [34]. The following equation modified from Wilkerson et al [14] was used to describe Model I: where PSA is the tumor burden represented by the PSA level at time t (in days); BAS is the baseline PSA; D is the drug effect on PSA, which decays over time due to the development of drug resistance; and G is the rate of growth of PSA due to proliferating tumor. T represents the time of PSA growth before the initiation of LHRH treatment; t = t if t ≤ time of the first LHRH dose (t) and t = t if t>t. T represents the time after the initiation of LHRH treatment; t = 0 if t ≤t and t = t-t if t>t.

It is widely recognized that tumor is highly heterogeneous, thus drug resistance can occur from treatment-induced selection of a subpopulation of drug-resistant cancer cells that existed in the original tumor [34, 35]. Therefore, Model II in Eqs 2 and 3 was used to describe this type of drug resistance: where PSA and PSA represent PSA levels at time t produced by drug-resistant and drug-sensitive cancer cell population, respectively. R is the proportion of drug-resistant cancer cells in the initial tumor population. R is expressed as e in the modeling process in order to obtain values of R between 0 and 1. G and G are growth rates of PSA due to replicating drug-resistant and drug-sensitive cancer cells, respectively. D is the drug effect on PSA levels due to killing of the drug-sensitive cancer cells. T represents the time of PSA growth before the initiation of LHRH treatment; t = t if t ≤ t and t = t if t>t. T represents the time after the initiation of LHRH treatment; t = 0 if t ≤t and t = t-t if t>t.

Inter-individual variability in the studied population was modeled for all parameters as follows: where θ is a model parameter, θ is the typical value of the corresponding parameter in the population, and η is a normally distributed random effect with a mean of 0 and a variance of ω. The residual error model in the analysis was additive error as follows: where DV and stand for j log observed and predicted concentration for the i subject. ε follows normal distribution with a mean of 0 and a variance of . Proportional error model and mixed error model were also assessed in the analysis.

Data were analyzed by nonlinear mixed effects modeling utilizing the Monte Carlo Parametric Expectation Maximization (MCPEM) method in NONMEM software (version 7.3, ICON Development). The details of the MCPEM algorithm have been presented elsewhere [36-38]. Briefly, two-stage hierarchical nonlinear mixed effects modeling was used to find the optimal population mean μ and variance Ω that best describe the observed data. Final population parameters μ and Ω were obtained by first evaluating the conditional mean and conditional variance for each subject using fixed values of μ and Ω (the expectation step E) according to Eqs 6 and 7, followed by evaluating updates to μ and Ω using Eqs 8 and 9 (the maximization step M) [37-39].

Expectation (E) step:

Maximization (M) step: where l(θ) is the likelihood function for subject i regarding to parameter θ given data; h(θ) is the density function of θ given μ and Ω. The weight W depends on likelihood function, l(θ), density function (suppressing dependence on μ and Ω), and the method of Monte-Carlo method used, for N randomly generated parameter vectors of θ. m represents the total number of subjects in the analysis. E and M steps were repeated until μ and Ω no longer change [36, 37]. At this point, final population parameters μ and Ω that best described the data were obtained. In MCPEM, the Monte-Carlo integration method was used to evaluate and during the expectation step [36, 37].

The parameters of the developed models were obtained by fitting the models to PSA levels from all subjects simultaneously. Comparison of alternative nested structural models was based on the typical goodness-of-fit diagnostic plots and likelihood ratio test [40, 41]. When comparing alternative nested hierarchical models, the differences in objective function value (OFV), defined as -2log(likelihood), are approximately chi-square distributed with n degrees of freedom (n is the difference in the number of parameters between the full and the reduced model) [41]. Differences in OFV of greater than 10.83 for one degree of freedom, corresponding to a significance level of 0.001, were used to discriminate two nested hierarchical models. This stringent criterion was used because of multiple comparisons inherent in the model selection procedure and random noise associated with the Monte-Carlo sampling technique employed in the MCPEM algorithm [39, 40, 42]. The number of simulated Monte-Carlo parameter sets (ISAMPLE) used for the E step evaluation determined the random noise associated with model parameters and likelihood estimation in the MCPEM method. Larger ISAMPLE was associated with lower random noise of likelihood for more reliable model selection but at the expense of increased computation times. Therefore, to achieve the desirable balance between random noise of model parameters/likelihood and computation time, a sequential approach with different ISAMPLE in different stages of the MCPEM was used in the analysis. First, ISAMPLE of 1000 for first 150 EM iterations and then ISAMPLE of 2000 for another 50 EM iterations were used in the burn-in phase of MCPEM to achieve the stationary phase rapidly. Then, ISAMPLE of 50,000 for 50 EM iterations was used to obtain stable and desirable final model parameter estimates using MCPEM model convergence criteria described previously [42]. In brief, last five EM iterations of population mean parameters and inter-subject variances with ISAMPLE of 50,000 at stationary state were examined using linear regression analysis. A Bonferroni method using the following equation was used to adjust the p value of linear regression analysis for multiple hypothesis testing: where α, α, and B were the Bonferroni-adjusted p value, preset p value (set at 0.05), and the number of tested parameters, respectively. If changes in all tested parameters across iterations were not statistically different from zero, then model convergence is assumed. To compare non-nested models such as Models I and II in this study, the following Akaike information criterion (AIC) was used to select the best model with lowest AIC value: where N is the total number of model parameters.

PSA concentrations below 0.1 ng/ml (lower limit of assay quantification) were treated as fixed point censored observations, and the maximum likelihood was used to fit the model to the censored observations [42, 43]. In this case, the likelihood for all data is maximized with respect to model parameters, and the likelihood for a censored concentration was taken to be the likelihood that the censored observation is truly below the limit of quantification. This approach allowed PSA levels <0.1 ng/mL to be included into model development to better characterize the PSA kinetics and treatment effect of leuprorelin.

Covariate analysis

In order to investigate and quantify relationships between important model parameters and patient-specific factors (covariates), final covariate model was developed using a modified Wald’s approximation method with backward elimination (WAM-BE) proposed by our group [44, 45]. In brief, all potential covariate-parameter relationships were incorporated into the best structural model to form the full model with covariates. Parameter estimates and the covariance matrix (COV) from the full model fit were used to calculate the Wald’s approximation statistics. Assuming the vector of fixed-effect covariate parameters was k×1 vector θ, and the corresponding COV for these parameters was k×k matrix C. θ could be partitioned to p×1 vector θ and q×1 vector θ, where θ were fixed-effect covariate parameters that need to be estimated and θ were covariate parameters restricted to zeros under the hypothesized submodel. Then θ and the corresponding C could be defined as follows:

The Wald’s approximation to the likelihood ratio test (LRT) statistic (Λ´) for the hypothesized submodel that θ = 0 was described by the equation:

It has been shown that maximum likelihood estimates of θ followed asymptotically multivariate normal distribution with covariance matrix C [46]. Under H0: θ = 0, the asymptotical distribution of Λ´ followed distribution, and hypothesis testing could be applied to decide whether or not to reject H0. The calculation of Λ´ in Eq 13 did not require running the submodels, and therefore, multiple hypotheses could be efficiently tested by using results from the full model fit. In WAM-BE method, the backward elimination (BE) process was used to screen and eliminate insignificant covariates from the full model based on the difference in values of Λ´between two models. Differences in Λ´ of greater than 10.83 for 1 degree of freedom, corresponding to a significance level of 0.001, were used to discriminate two nested hierarchical models. The best models selected by BE approach based on Λ´ values were then used as the starting models for selecting the final covariate model using actual NONMEM runs with the BE process.

Continuous covariates were added into the structural model as follows [40]: where COV is the value of j continuous covariate, and θ is the effect the j covariate has on the i parameter.

Categorical covariates were added as follows: where COV is an indicator variable (e.g., if the categorical variable is gender, COV = 1 when the subject is male, and COV = 0 when the subject is female).

The final model was evaluated by routine diagnostic plots. Additionally, prediction-corrected visual predicted check (VPC) was done in PsN (ver. 4.9.0, Uppsala University, Sweden) using 1000 simulations and automatic binning, to assess the performance of the final model [47]. The final population model was used to examine the effects of statistically significant covariates on clinical outcomes associated with leuprorelin treatment. PSA progression is defined by the Prostate Cancer Clinical Trials Working Group (PCWG2 and PCWG3) as “the date that a 25% or greater increase and an absolute increase of 2 ng/mL or more from the nadir is documented, which is confirmed by a second value obtained 3 or more weeks later” [48, 49]. The percentages of leuprorelin-treated subjects with PSA progression within one, two, and three years simulated from 1000 subjects using 5th percentile, median and 95th percentile values of continuous covariates and different categories of categorical variables were obtained and compared. Simulation was conducted with the final population model using R.

The parallel computing platform of NONMEM was implemented in a single workstation equipped with dual Intel Xeon E5-2698 v4 20-cores CPU with 2.20 GHz, a Windows 10 Enterprise operating system and Intel Parallel Studio XE 2016 Fortran compiler.

Results

A total of 1113 PSA observations from 264 subjects were used for model development. Baseline demographic data and laboratory values of all subjects are summarized in Table 1. The study population had a median age of 80 and median baseline PSA of 8.5 ng/mL. There was no missing baseline demographic data and laboratory values among subjects included in the analysis.

Table 1

Summary of baseline demographic data and laboratory values of subjects included in model development.

Baseline characteristics	Median (range) or counts
Age (years)	80 (60–100)
Race
Caucasian	189
Black	59
Hispanic/other	16
Region
South	196
West	20
Midwest	42
Northeast	6
Antiandrogen use
Yes	33
No	231
aAST (IU/L)	20 (9–91)
bALT (IU/L)	18 (4–110)
cSCR (mg/dL)	1.10 (0.700–9.30)
dALP (IU/L)	76.5 (23.0–3640)
eALB (g/dL)	4.13 (2.90–4.80)
fHGB (g/dL)	13.6 (6.80–17.4)
Baseline gPSA (ng/mL)	8.50 (0.200–782)

Data are presented as median (range) or counts.

aaspartate transaminase

balanine transaminase

cserum creatinine

dalkaline phosphatase

ealbumin

fhemoglobin

gprostate-specific antigen.

Model I that assumed drug resistance was developed from PSA-producing PCa cells that were initially sensitive to LHRH treatment, had an AIC value of 865.432. Model II that assumed drug resistance was originated from treatment-induced selection of a subpopulation of drug-resistant PSA-producing cancer cells, had a lower AIC of 295.671. Therefore, Model II was selected as the structural model. Use of proportional or mixed error models resulted in unstable models, so an additive error model was used. HGB was a significant covariate that affects the proportion of drug-resistant cells (R or e) in the original tumor. Baseline PSA (BAS) and antiandrogen use (AND) were significant covariates on the drug killing effect of drug-sensitive PSA-producing cancer cell population (D). Addition of these three covariates into the model resulted in significant decrease in OFV (50.008; p<0.001 for degrees of freedom = 3). The final model with covariates was presented as follows: where θ and θ are model parameters used to describe the effect of HGB on RP and BAS on D, respectively. θ is used to describe the effect of AND on D. IND = 1 if the patient used antiandrogen within 30 days of leuprorelin initiation and = 0 otherwise. The parameter estimates of the final model with covariates are presented in Table 2. All model parameters were estimated with good precision (percent coefficient of variation, %CV<50). The population estimate of D was 3.78 x 10−2 day-1. Population estimate of growth rate for drug-sensitive PSA-producing tumor cell population (G) was 1.96 x 10−3 day-1, corresponding to a PSA doubling time of 354 days. Population estimate of RP was 3.94, indicating that 1.94% (calculated from R = e) of the original PSA-producing cancer cell population was inherently resistant to leuprorelin treatment. Population estimate of growth rate for drug-resistant PSA-producing cancer cell population (G) was 6.54 x 10−4 day-1, corresponding to a PSA doubling time of 1060 days.

Table 2

Parameter estimates of the final PSA kinetics model.

Parameter	Estimate	a%CV
Structural Model
bDS (day-1)	3.78 x 10−2	6.19
cGs (day-1)	1.96 x 10−3	22.5
dRP	3.94	7.44
eGR (day-1)	6.54 x 10−4	28.4
Interindividual Variability
ωDS2	0.453	14.8
ωGS2	2.59	19.7
ωRP2	0.944	16.4
ωDR2	3.76	21.1
Covariate Model
fHGB on RP (θHGB_RP)	2.30	24.7
gBAS on DS (θBAS_DS)	0.174	24.5
hAND on DS (θAND_DS)	0.677	30.0
Residual Variability
Additive error (σadd)	2.01x10-1	3.49

a%CV = percent coefficient of variation

bdrug effect on drug-sensitive tumor cells

cgrowth rate of drug-sensitive cells

dexp(-RP) represents the fraction of drug-resistant tumor cells in the original tumor

egrowth rate of drug-resistant cells

feffect of hemoglobin level on RP

geffect of baseline prostate-specific antigen level on Ds

heffect of antiandrogen use on Ds

Large between-patient variability was observed in these parameters, possibly due to diverse PSA kinetics in the studied population. The estimated value of θ was 2.30, suggesting that subjects with lower HGB levels had lower RP and therefore were less likely to respond to leuprorelin treatment. Typical values of RP were 2.37, 3.94, and 5.73 for subjects with 5th percentile (10.9 g/dL), median (13.6 g/dL), and 95th percentile (16.0 g/dL) levels of HGB, respectively. The corresponding proportions of resistant PSA-producing cancer cell population in the original tumor were 9.36%, 1.94%, and 0.326%, respectively. The estimated value of θ was 0.174, indicating that higher baseline PSA levels were associated with better LHRH-treatment effect on PSA level due to killing of drug-sensitive cancer cell population. Typical values of D were 2.32 x 10−2, 3.78 x 10−2, and 5.99 x 10−2 day-1 for subjects with 5th percentile (0.515 ng/mL), median (8.50 ng/mL), and 95th percentile (120 ng/mL) levels of BAS, respectively. The estimated value of θ was 6.77 x 10−1, implying that subjects with antiandrogen use within 30 days of leuprorelin initiation was associated with 96.8% higher LHRH treatment effect on killing of drug-sensitive PSA-producing cancer cell population compared to those without antiandrogen use.

As shown in Fig 1, diverse PSA kinetics profiles from different patients were well described by the final model. The diagnostic plots in Fig 2 and the prediction-corrected VPC plot shown in S2 Fig show that the observed PSA values were generally in agreement with predicted PSA levels, and that the final model provided a good fit for the data. The final model was then used in simulation to assess effects of HGB, BAS and AND on PSA progression within one, two and three years, and the results are shown in Fig 3. The percentages of leuprorelin-treated subjects with PSA progression within one year for subjects with 5th percentile, median and 95th percentile of HGB were 19.8, 13.9 and 10.5, respectively (Fig 3A). 38.8, 28.2, and 22.1% of the leuprorelin-treated subjects with 5th percentile, median and 95th percentile of HGB were expected to experience PSA progression within three years after initiation of leuprorelin treatment. The percentages of leuprorelin-treated subjects with PSA progression within one year in subjects with 5th percentile, median and 95th percentile of BAS were 6.9, 13.9, and 25.6, respectively (Fig 3B). The percentages of leuprorelin-treated subjects with PSA progression within three years in subjects with 5th percentile, median and 95th percentile of BAS were 15.8, 28.2, and 47.8, respectively (Fig 3B). 28.5% of leuprorelin-treated subjects with antiandrogen use and 28.2% of those without antiandrogen use were expected to experience PSA progression within three years of treatment period.

Fig 1

Profiles of PSA kinetics in representative patients with hormone-sensitive prostate cancer.

Black line represents model-predicted PSA levels. Open circles and solid circles represent observed PSA levels above the lower limit of quantification (LLOQ) and observed PSA levels below the LLOQ, respectively. Red and green vertical dashed lines represent the first and the last recorded date of leuprorelin treatment, respectively. Blue horizontal dotted line represents the LLOQ value of the PSA assay. Black arrows represent the fill dates of leuprolide.

Fig 2

Diagnostic plots for the final disease progression model.

From left to right and top to bottom: Observed log PSA concentrations (DV) versus individual log predicted values (IPRED), individual weighted residuals (IWRES) versus the IPRED, IWRES versus time, and expected conditional weighted residuals (ECWRES) versus the expected log predicted values (EPRED). The red line represents the loess regression line.

Fig 3

Simulation results for percentages of subjects with PSA progression within one, two, and three years.

3A and 3B show simulated PSA progression in subjects with 5th percentile, median and 95th percentile of hemoglobin level and baseline PSA level, respectively.

Discussion

To the best of our knowledge, this was the first study to use medical claims data to construct a population-based disease progression model to describe PSA kinetics in patients with hormone-sensitive PCa treated with leuprorelin. Furthermore, the model could quantify effects of patient-specific factors (covariates) on model parameters and PSA kinetics in these patients. Although leuprorelin was selected as the only form of ADT in this analysis, we would expect that our results could be extrapolated to other LHRH agonists, given the similarity in their efficacy reported in clinical trials [50, 51].

Population data analysis using nonlinear mixed effects modeling approach has been shown to provide precise and robust estimates of drug-related model parameters and their population variability using clinical data [41, 52]. Among different nonlinear mixed effects modeling methods, MCPEM was used to successfully develop the final population disease progression model in this study, due to the ability of this method to precisely approximate the true log-likelihood, and its stability and amenability to efficient parallel computing when dealing with complex models with large datasets [38, 39, 53].

One of the important goals in population-based disease progression modeling is to quantify effects of covariates on model parameters and clinical outcomes. Established quantitative relationships are useful for helping clinicians make conscious decisions on selecting optimal drug treatment, with the ultimate goal of implementing individualized or personalized medicine [54]. However, development of population-based disease progression models with covariate-parameter relationships is both time-consuming and labor-intensive, as many submodels with different covariate-parameter relationships need to be created, tested, and compared. There is a total of 24x11 = 1.76 x 1013 possible submodels with different covariate-parameter relationships in this study. Due to the structural model complexity and large number of observations, a single covariate model run with ISAMPLE of 50,000 executed in parallel mode with 25 Intel Xeon 2.2 GHz E5-2968 v4 computing cores took about 10 min to complete. Therefore, it was impossible to explore all the possible submodels in this setting. Hence, the WAM-BE method originated by our group was used to develop covariate models [45]. WAM-BE only requires a single full model fit and does not require fitting the submodels to estimate the difference in OFV between two tested models for covariate selection. In addition, WAM-BE is designed to overcome the inherent limitation of the original WAM method in population data analysis [44] and it achieved comparable results with significantly shorter computational time in models with large numbers of model parameters and tested covariates. Although the original WAM method uses the Wald approximation to LRT statistic (Λ´) from the full model instead of actual NONMEM model runs to screen covariate models, it is extremely inefficient in the presence of large numbers of tested model parameters and covariates due to the need to calculate Λ´ values for all possible covariate models. In this study, 1.76 x 1013 possible covariate models needed to be screened and tested with the original WAM method. Assuming the time to compute Λ´ value for a covariate model was about 10−6 seconds, a total of 1.76 x 107 seconds or 203 days would be needed to calculate Λ´ values of all possible covariate models. To overcome this limitation, we used the sequential backward elimination (BE) process to efficiently eliminate any insignificant covariate models to generate the best starting models, which were subject to subsequent actual NONMEM runs with the BE process to select the final covariate model. Therefore, the WAM-BE method significantly reduced the number of covariate models needed to be tested using WAM-derived Λ´ values and actual NONMEM runs, and it only took about 39 min for our workstation to develop the final covariate model.

ADT is the standard of care for patients with hormone-sensitive PCa, but most patients develop castration resistance within 1–3 years [53]. The underlying mechanism is thought to be multifactorial and involve several molecular and genetic alterations, and the “adaptation” and the “clonal selection” models have been proposed to explain this phenomenon [55]. With the “adaptation” model, early PCa cells are assumed to have similar androgen requirement for survival and growth, and castration resistance stems from genetic or epigenetic alterations in the cells. With the “clonal selection” model, it is assumed that early PCa cells have heterogeneous androgen requirement, and after ADT castration-resistant cells gradually outgrow castration-sensitive cells due to the survival advantage of the former [56].

In our model, PSA was used as a marker of tumor burden in the modeling of disease progression, as PSA levels are routinely used for disease monitoring and surveillance purposes in PCa [20-23], and they are readily retrievable from health claims database. PSA level has been shown to be closely related to the tumor volume [57], and it was found to reflect the androgen milieu in patients with localized PCa on ADT [51]. A rising PSA is usually the first sign of tumor regrowth, followed by worsening of disease identified by imaging and development of clinical symptoms [49]. In this study, Model I was the mathematical representation of “adaptation” model that assumed drug resistance was developed from PSA-producing cancer cells initially sensitive to LHRH treatment. On the other hand, Model II was the “clonal selection” model that assumed drug resistance was developed from a subpopulation of drug-resistant PSA-producing cancer cells in the original tumor. Our predefined model selection criteria determined that Model II was superior than Model I at describing the PSA profiles of our study subjects, and this finding was consistent with the observations from some recent studies [55].

Based on recent recommendations from the PCWG, PSA responses to drug treatment that affect tumor cell kill can be categorized into three distinct patterns: 1) significant and sustained PSA decline after drug treatment for “responders”, 2) initial decline followed by a slow PSA increase for “partial responders”, and 3) a transient decrease followed by a rapid PSA increase observed in “non-responders” [49]. As shown in Fig 1, the final model in this study was able to adequately describe these three distinct PSA responses to leuprorelin treatment. Furthermore, the model was able to provide additional biological insight into the observed PSA patterns. For example, Subject C in Fig 1C showed least response to treatment and most significant increase in PSA within the follow-up period among the representative subjects. The modeling results suggested that this subject had a high percentage of drug-resistant cancer cell population of 28.1% and high growth rate of the resistant tumor cell population of 7.97 x 10−3 day-1, both of which were more than 10 times higher than the population mean values. On the other hand, the profile of Subject B in Fig 1B was characterized by a decrease followed by a slow increase of PSA. As shown in Fig 4B, the rapid but transient decline in the total PSA level after leuprorelin treatment was due to the killing of drug-sensitive PSA-producing cancer cell population. However, the high growth rate of drug-resistant PSA-producing cancer cell population of 1.53 x 10−2 day-1 (compared to the population estimate of 6.54 x 10−4 day-1) contributed to the subsequent steady rise in PSA level observed in this subject. On the other hand, PSA profile of Subject A in Fig 1A was characterized by significant and sustained PSA decline after leuprorelin treatment and therefore could be described as a "responder”. The low percentage of drug-resistant PSA-producing cancer cell population of 6.36 x 10−2% combined with the slow growth rate of the drug-resistant PSA-producing cancer cell population of 3.53 x 10−4 day-1 contributed to the sustained PSA decline in this subject (Fig 4A).

Fig 4

Simulated total PSA and PSA produced by drug-sensitive and drug-resistant tumor cells in representative subjects.

4A and 4B demonstrate simulated PSA levels in a “responder” Subject A and a “partial responder” Subject B, respectively. Black line represents model predicted total PSA level. Pink and grey solid lines represent model predicted PSA from drug-sensitive and drug-resistant tumor cells, respectively. Red and green vertical dashed lines represent the first and the last recorded date of leuprorelin treatment, respectively. Blue horizontal dotted line represents the lower limit of quantification (LLOQ) of the PSA assay.

Only a few studies have attempted to develop predictor models for development of castration resistance in patients with hormone-sensitive PCa. Most of these predictive models were developed using multivariate analysis that provides very little biological insight on the relationship between identified predictors and disease progression [8-12]. Several PSA kinetic parameters such as time to PSA nadir and PSA nadir level have been shown to be significant predictors of disease progression for PCa patients receiving ADT [9-11]. However, the complex quantitative relationship between clinical predictors and PSA kinetics and disease progression in ADT-treated PCa patients remains largely unknown. In this study, we developed the first mechanistic population-based disease progression model that allowed us to investigate such complex relationship. Baseline PSA and HGB were identified as important covariates on PSA kinetics, consistent with previous findings that these two covariates played an important role in the development of castration resistance [8, 9, 12]. However, these two covariates affected different model parameters and hence had very different results on PSA kinetics and biochemical progression. Both baseline PSA and antiandrogen use were statistically significant covariates on drug-killing effect on leuprorelin-sensitive PSA-producing cancer cells, and HGB significantly affected the fraction of cancer cells resistant to leuprorelin treatment. Drug-killing rates on leuprorelin-sensitive PSA-producing cancer cells were higher in subjects with higher baseline PSA levels and those with antiandrogen use within the 30 days of leuprorelin treatment initiation, compared to those with lower baseline PSA levels and those without antiandrogen use. Subjects with lower hemoglobin levels had higher fraction of leuprorelin-resistant cancer cells in the original tumor.

Subsequently, simulation was performed to evaluate the relationship between these covariates and PSA kinetics and disease progression. The results showed that subjects with higher baseline PSA levels were more likely to experience PSA progression (Fig 3B), which may seem counterintuitive given that higher baseline PSA was associated with higher drug-killing effect on leuprorelin-sensitive PSA-producing cancer cells in our model. Further analysis of the simulated results showed that the baseline PSA level had a complicated effect on the overall PSA kinetics. The simulated median time to nadir levels was 294, 198, and 135 days after initiation of leuprorelin in subjects with 5th percentile (0.515 ng/mL), median (8.50 ng/mL) and 95th percentile (120 ng/mL) of baseline PSA levels, respectively. The simulated median nadir levels were 0.0220, 0.290, and 3.47 ng/mL for leuprorelin-treated subjects with 5th percentile, median and 95th percentile of baseline PSA levels, respectively. Therefore, subjects with higher baseline PSA levels had higher velocity of PSA decline and shorter time to achieve the nadir level after treatment initiation, which was consistent with higher drug-killing effect on leuprorelin-sensitive cancer cells in these subjects. Furthermore, subjects with higher baseline PSA levels also had higher nadir PSA levels, possibly due to presence of higher numbers of leuprorelin-resistant cancer cells, and therefore they were more likely to experience PSA progression. Overall, our results confirmed findings from previous studies that associated higher baseline PSA levels with shorter time to castration resistance, and they were also in agreement with observations from Ji et al that higher PSA nadir, higher velocity of PSA decline and shorter time to PSA nadir were predictive of increased risk of progression to castration resistance [11]. On the other hand, subjects with lower levels of HGB were shown to have a higher chance of PSA progression (Fig 3A), consistent with results from previous studies [8, 12]. Mechanistically, the simulated median time to nadir levels were 155, 198, and 240 days after initiation of leuprorelin in subjects with 5th percentile (10.9 g/dL), median (13.6 g/dL) and 95th percentile (16.0 g/dL) of HGB levels, respectively. The simulated median nadir levels were 1.19, 0.290, and 0.0574 ng/mL for subjects with 5th percentile, median and 95th percentile of HGB levels, respectively. Again, we could see that subjects with lower HGB levels had higher PSA nadir, higher velocity of PSA decline and shorter time to PSA nadir, which were the identified PSA kinetic risk factors predictive of shorter time to biochemical progression. Lastly, while use of antiandrogen within 30 days of initiation of leuprorelin treatment was identified as a statistically significant covariate affecting drug-killing effect on leuprorelin-sensitive PSA producing tumor cells, it had minimum effect on PSA progression in our simulation. One plausible explanation was that as antiandrogen use could only be temporary in the study population, the observed inhibitory effect was most likely on PSA flare alone, which has not been found to be associated with tumor progression and negative outcomes [58].

There are several limitations to our study. First, while medical claims database can be used for longitudinal population studies to address questions in real-life clinical practice, the highly heterogeneous nature of the database posed technical challenges to data analysis. For example, patients on both continuous and intermittent regimens of leuprorelin were included in the analysis, as it was difficult to quantitatively identify and separate these two types of treatment regimens based on claims data alone. Additionally, because of the heterogeneity of dosing regimens in the studied population, patient compliance could not be accurately assessed. Intermittent ADT has been extensively tested in PCa patients since the 1980s, mainly due to evidence of its anticancer efficacy and its ability to reduce ADT-related adverse effects [59]. With intermittent ADT, the patient receives ADT for several months consecutively until PSA levels drop below a predetermined threshold, after which the patient enters an off-treatment period [60]. Patients receiving the two treatment modalities may have different patterns of disease progression. Nevertheless, two randomized trials showed that intermittent ADT was noninferior to continuous ADT in terms of survival endpoints in PCa patients [61, 62], while the largest trial so far comparing continuous and intermittent ADT could not conclude or exclude noninferiority of intermittent ADT [63]. In our preliminary analysis, dose intensity (leuprolide dose received normalized to an expected dose of 7.5 mg per month) was calculated. As expected, dose intensity varied greatly among patients on continuous and intermittent ADT, and it did not have statistically significant effect on model parameters. Additionally, exclusion of a large proportion of patients (86.4%) diagnosed with malignant PCa and treated by leuprorelin in our effort to build an explainable model, as well as the presence of large between-patient variability observed in some model parameter estimates, may pose a limitation to interpretation of our findings. Nevertheless, retrospective nature of this study suggested that our findings will need to be validated in prospective clinical studies. Second, some clinical measures that assess disease severity and/or metastasis spread, including metastasis stage, the Gleason score, Eastern Cooperative Oncology Group (ECOG) performance status and the Soloway score, were shown to be predictive of biomedical progression [9, 11, 12]. However, they were not available in the medical claims database, and therefore they could not be evaluated for their importance in our model. Similarly, though LHRH agonists have been assumed to cause testosterone suppression [64], testosterone levels were not available in our patient database for us to investigate the possibility of incomplete testosterone suppression and its potential relationship with PSA kinetics. Third, as survival data were not available in our database, analyses linking PSA kinetics and biomedical progression to survival endpoints were not conducted in this analysis. Last but not least, though PSA level has been previously shown to be significantly associated with tumor volume [57], it is not as accurate as quantitative imaging biomarkers in representing the tumor burden, and the PSA response does not always associate with treatment response or survival [65].

Conclusions

In summary, in this study a novel population-based modeling approach was used to provide the first mechanistic insight on how PSA kinetics contributed to the subsequent development of castration resistance in patients with hormone-sensitive PCa treated with leuprorelin. Compared to previous studies, our study successfully modeled the underlying resistance mechanism of PCa cells and provided new biological insight into PSA kinetics in this patient population of interest. The application of population-based disease progression model to existing clinical data allowed estimation of tumor resistant patterns and growth/regression rates that could greatly enhance our understanding of how hormone-sensitive PCa responds to LHRH agonists. It may serve as a platform for incorporating more comprehensive health data in the future, including laboratory measurements, genomic and proteomic data to further personalized medicine in patients with hormone-sensitive PCa.

Patient selection flow diagram.

(TIF)

Click here for additional data file.

Prediction-corrected visual predictive check of the final model.

Blue circles represent prediction-corrected observations (log-transformed PSA concentrations). The solid red line represents the median prediction-corrected observations, and the semitransparent red area represents the simulated 95% confidence interval for the median. The dashed red lines represent the 5% and 95% percentiles of prediction-corrected observations, and the semitransparent blue areas represent their respective simulated 95% confidence intervals.

(TIF)

Click here for additional data file.

Data dictionary for the Humana dataset.

(XLSX)

Click here for additional data file.

4 Dec 2019

PONE-D-19-27455

Using medical claims database to develop a population disease progression model for leuprorelin-treated subjects with hormone-sensitive prostate cancer

PLOS ONE

Dear Dr. Ng,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

==============================

Please follow the suggestions of the two reviewers who are experts in prostate cancer field and revise the manuscript to answer their questions.

==============================

We would appreciate receiving your revised manuscript by Jan 18 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Chih-Pin Chuu, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you have indicated that data from this study are available upon request. PLOS only allows data to be available upon request if there are legal or ethical restrictions on sharing data publicly. For information on unacceptable data access restrictions, please see http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions.

In your revised cover letter, please address the following prompts:

a) If there are ethical or legal restrictions on sharing a de-identified data set, please explain them in detail (e.g., data contain potentially identifying or sensitive patient information) and who has imposed them (e.g., an ethics committee). Please also provide contact information for a data access committee, ethics committee, or other institutional body to which data requests may be sent.

b) If there are no restrictions, please upload the minimal anonymized data set necessary to replicate your study findings as either Supporting Information files or to a stable, public repository and provide us with the relevant URLs, DOIs, or accession numbers. Please see http://www.bmj.com/content/340/bmj.c181.long for guidelines on how to de-identify and prepare clinical data for publication. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories.

We will update your Data Availability statement on your behalf to reflect the information you provide.

3. In the ethics statement in the manuscript and in the online submission form, please provide additional information about the patient records used in your retrospective study, including: the date range (month and year) during which patients' medical records were accessed.

4.  Thank you for stating the following in the Financial Disclosure section:

The authors received no specific funding for this work.

We note that one or more of the authors are employed by a commercial company: NewGround Pharmaceutical Consulting LLC

1. Please provide an amended Funding Statement declaring this commercial affiliation, as well as a statement regarding the Role of Funders in your study. If the funding organization did not play a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript and only provided financial support in the form of authors' salaries and/or research materials, please review your statements relating to the author contributions, and ensure you have specifically and accurately indicated the role(s) that these authors had in your study. You can update author roles in the Author Contributions section of the online submission form.

Please also include the following statement within your amended Funding Statement.

“The funder provided support in the form of salaries for authors [insert relevant initials], but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the ‘author contributions’ section.”

If your commercial affiliation did play a role in your study, please state and explain this role within your updated Funding Statement.

2. Please also provide an updated Competing Interests Statement declaring this commercial affiliation along with any other relevant declarations relating to employment, consultancy, patents, products in development, or marketed products, etc.

Within your Competing Interests Statement, please confirm that this commercial affiliation does not alter your adherence to all PLOS ONE policies on sharing data and materials by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests) . If this adherence statement is not accurate and  there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

Please include both an updated Funding Statement and Competing Interests Statement in your cover letter. We will change the online submission form on your behalf.

Please know it is PLOS ONE policy for corresponding authors to declare, on behalf of all authors, all potential competing interests for the purposes of transparency. PLOS defines a competing interest as anything that interferes with, or could reasonably be perceived as interfering with, the full and objective presentation, peer review, editorial decision-making, or publication of research or non-research articles submitted to one of the journals. Competing interests can be financial or non-financial, professional, or personal. Competing interests can arise in relationship to an organization or another person. Please follow this link to our website for more details on competing interests: http://journals.plos.org/plosone/s/competing-interests

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This study well demonstrated pharmacokinetic in leuprorelin treatment of prostate cancer from insurance population data.

The authors are trying to mechanistically explain leuprorelin's effect on prostate cancer cell proliferation by PSA kinetics, which may be quite ambitious

But the article did find itself a role in explaining tumor resistant patterns, which can be of great value in selected patients to early modify prostate cancer treatment

Advantages

1. Proposed a good model of drug resistance development by demonstrating individual leuprorelin sensitive and resistant cancer cell growth rate using PSA findings.

2. Reconfirmed important covariates in prostate cancer development including baseline PSA and hemoglobin, clarify important factors associated with leuprorelin treatment

3. Specify the portion of leuprorelin-resistant cancer cell population and its growth rate can improve patient’s follow-up and the early detection of most at risk patients.

4. The data source is unique and stands for the real-world population. No pharmacokinetics of leuprorelin researches with Humana database to date yet.

5. Provide population mean prostate growth rate and resistance portion for future reference.

Shortcomings

1. Large between-patient variability may indicate that few models may not best estimate leuprorelin on drug resistance of prostate cancer,

2. Need to exclude majority of patients (2035/2354 =86.4%) under leuprorelin treatment in order to build up a explanable model, indicated that the model may not be universally applicable.

3. Figure legends can be refined for better understanding eg. Fig.3A legends should add HGB and Fig.3B BAS on the figure.

4. This research did not take into account the incomplete testosterone suppression and its effect on PSA kinetics.

Reviewer #2: This manuscript is a study which looks at disease progression model for patients with hormone sensitive prostate cancer on leuprorelin. The conclusion is that the proportion of the original prostate cancer cells inherently resistant to treatment was estimated to be 1.94%. This is an interesting study and it can be published after some revision.

Strengths of this study:

1 Large sample size.

2 Well data analyzed.

This paper is easy to follow and well structured. But the characters of the prostate cancer patients should be clarified. The prostate cancer staging, Gleason score of prostate cancer and comorbility data are not available in the paper. This is a major limitation of the study. The authors should tell us the distribution of prostate cancer risk factors in these patients.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

17 Feb 2020

Response to Reviewers

Reviewer #1 Comments

This study well demonstrated pharmacokinetic in leuprorelin treatment of prostate cancer from insurance population data

The authors are trying to mechanistically explain leuprorelin's effect on prostate cancer cell proliferation by PSA kinetics, which may be quite ambitious

But the article did find itself a role in explaining tumor resistant patterns, which can be of great value in selected patients to early modify prostate cancer treatment

We would like to thank the reviewer’s comments, and we have made some changes to our manuscript accordingly.

Advantages

1. Proposed a good model of drug resistance development by demonstrating individual leuprorelin sensitive and resistant cancer cell growth rate using PSA findings.

2. Reconfirmed important covariates in prostate cancer development including baseline PSA and hemoglobin, clarify important factors associated with leuprorelin treatment

3. Specify the portion of leuprorelin-resistant cancer cell population and its growth rate can improve patient’s follow-up and the early detection of most at risk patients.

4. The data source is unique and stands for the real-world population. No pharmacokinetics of leuprorelin researches with Humana database to date yet.

5. Provide population mean prostate growth rate and resistance portion for future reference.

Shortcomings

1. Large between-patient variability may indicate that few models may not best estimate leuprorelin on drug resistance of prostate cancer

We acknowledge that the presence of large between-patient variability is a limitation of our study, and we have included this in our discussion section L621-625.

2. Need to exclude majority of patients (2035/2354 =86.4%) under leuprorelin treatment in order to build up a explanable model, indicated that the model may not be universally applicable.

We acknowledge that the need to exclude the majority of patients is a limitation of our study, and we have included this our discussion section L621-625, which reads, “Additionally, exclusion of a large proportion of patients (86.4%) diagnosed with malignant PCa and treated by leuprorelin in an effort to build an explainable model, as well as the presence of large between-patient variability observed in some model parameter estimates, may pose a limitation to interpretation of our findings.”

3. Figure legends can be refined for better understanding eg. Fig.3A legends should add HGB and Fig.3B BAS on the figure.

Figures 3A and 3B have been corrected to include HGB and BAS on the figures.

4. This research did not take into account the incomplete testosterone suppression and its effect on PSA kinetics.

Testosterone levels are not available in our dataset, as it was not routinely monitored during the period covered by our database (2007-2011). We have acknowledged this limitation in our discussion section L631-634, which reads, “Similarly, though LHRH agonists have been assumed to cause testosterone suppression [63], testosterone levels were not available in our patient database for us to investigate the possibility of incomplete testosterone suppression and its relationship with PSA kinetics.”

Reviewer #2: This manuscript is a study which looks at disease progression model for patients with hormone sensitive prostate cancer on leuprorelin. The conclusion is that the proportion of the original prostate cancer cells inherently resistant to treatment was estimated to be 1.94%. This is an interesting study and it can be published after some revision.

We would like to thank the reviewer’s comments, and please see our response below.

Strengths of this study:

1 Large sample size.

2 Well data analyzed.

This paper is easy to follow and well structured. But the characters of the prostate cancer patients should be clarified. The prostate cancer staging, Gleason score of prostate cancer and comorbility data are not available in the paper. This is a major limitation of the study. The authors should tell us the distribution of prostate cancer risk factors in these patients.

The prostate cancer staging and Gleason score of the patients are, unfortunately, not available in our dataset, as they were not associated with any diagnostic/procedural/laboratory codes and therefore were not captured by de-identified data extraction from the electronic health records. These limitations were acknowledged in the discussion section L626-631.

Though comorbidities influence survival in patients with metastatic prostate cancer because they are associated with non-cancer associated deaths, to our knowledge, the effect of comorbidities is on biomedical progression (i.e., progression to castration resistance) has not been reported in the literature. As we focused on biomedical progression in castration-sensitive prostate cancer patients instead of survival in patients with castration-resistant disease, we think that comorbidities are likely not as important than the other factors mentioned and studied in our manuscript, including demographics and laboratory values. Additionally, electronic health record (EHR) problem lists have been shown to have poor sensitivity for detecting major comorbidities compared to free-text notes in the EHR (Daskivich TJ et al, Am J Manag Care. 2018 Jan 1;24(1):e24-e29); however, the latter was not accessible to the authors due to the de-identified nature of the data source.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

4 Mar 2020

Using medical claims database to develop a population disease progression model for leuprorelin-treated subjects with hormone-sensitive prostate cancer

PONE-D-19-27455R1

Dear Dr. Ng,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Chih-Pin Chuu, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

9 Mar 2020

PONE-D-19-27455R1

Using medical claims database to develop a population disease progression model for leuprorelin-treated subjects with hormone-sensitive prostate cancer

Dear Dr. Ng:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Prof. Chih-Pin Chuu

Academic Editor

PLOS ONE