We systematically searched PubMed, EMBASE, and the Cochrane Library from inception to February 6, 2020, for randomized clinical trials (RCTs). For the database search, we used “ovar* and cancer*, carcinoma*, neoplasm, or tumour and Olaparib, LYNPARZA, rucaparib, RUBRACA, niraparib, ZEJULA, AZD2281, AZD-2281, KU-0059436, CO-338, AG14699, MK4827” as search terms in all fields. The search was restricted to articles published in English.

The inclusion criteria for the analysis were as follows: 1) studies on patients with ROC for maintenance treatment; 2) studies on drugs, including olaparib, rucaparib, and niraparib, approved by the FDA; 3) randomized studies and published progression-free survival (PFS) results of the PARPi arm; and 4) information on biomarkers (BRCA gene or HRD status).

The exclusion criteria were: 1) reviews, case reports, or meta-analyses; 2) conference abstracts; and 3) studies on non-maintenance therapy.

The following information was extracted using Microsoft Excel software (version 2016): 1) characteristics of the literature (authors, year of publication, DOI/PMID, and clinical trial registration number); 2) trial design (randomized information and blind information); 3) characteristics of participants (disease information, biomarker status, Eastern Cooperative Oncology Group (ECOG) status, and sample size); 4) treatment information including PARPi and combination drugs (generic drug name, dosage, dose unit, strength, and administration information); 5) PFS results.

All data were extracted from the included studies by two independent researchers. Any dispute was resolved through discussion with a third independent researcher. Digitization software Get Data (version 2.26) was used to obtain graphical data. When the extraction errors of the graphics data were higher than 2%, data extraction was repeated, and the mean values were used as the final results.

The quality of the literature was assessed based on the Cochrane risk-of-bias criteria by two independent investigators, and each item was graded as low, high, or unclear. In total, seven items were used to evaluate bias in each trial, including random sequence generation, allocation concealment, blinding of participants and personnel, blinding of outcome assessment, incomplete outcome data, selective reporting, and other biases. The term “other bias” included trials sponsored by drug companies or studies in which the baseline characteristics of the intervention groups were not similar.

The quality of the literature was graded independently by two researchers, and inconsistencies were resolved by a third independent researcher.

The primary endpoint was PFS, which was analyzed using the parametric survival function. Four basic hazard ratio models, including the exponential, Gompertz, Weibull, and log normal models (Eqs 1–4) (Holford, 2013; Ding et al., 2020), were evaluated. Once the hazard ratio model was defined, the PFS data were fitted using Eq. 5. The inter-study variability and residual error were also considered during the establishment of the basic model (Eqs 6–8). h ( t ) = λ 0 (1) h ( t ) = λ 0 × e β × t (2) h ( t ) = λ 0 × e β × ln ( t ) (3) h ( t ) = ( σ t 2 π ) − 1 e ( − 1 / 2 Z 2 ) 1 − φ ( z ) ,   Z = ln ( t ) − μ σ (4) S ( t ) = exp ( − ∫ 0 t h ( t ) d t ) (5) P i = P p o p × e η i (6) O b s j , i = P r e d j , i × ( 1 + S E j , i × ε ( 1 ) j , i ) + S E j , i × ε ( 2 ) j , i (7) S E j , i = O b s j , i × ( 1 − O b s j , i ) N i (8)

The hazard function h_(t) (Eqs 1–4) is the instantaneous hazard of dying at time t, and in Eqs 1–3 λ₀ and β represent the hazard rate at time 0 and change in the regression coefficient of the hazard rate over time, respectively. In Eq. 4, μ and σ are the median and standard deviation of the distribution, respectively. In Eq. 6, P_i is the individual prediction of the model parameter, and P_pop is the population prediction of the corresponding model parameter. η_i represents the inter-study variability, which is assumed to be normally distributed, with a mean of 0 and a variance of ω². In Eq. 7, Obs_j,i represents the observed PFS rate at time point j in study i. Pred_j,_i is the corresponding individual prediction of the final model. ɛ_(1)j,_i and ɛ_(2)j,_i are the proportional residual and added residual errors, respectively, at the time point j of study i, assumed to be normally distributed with a mean of 0 and variance of σ₁ ² and σ₂ ². The residual variability was weighed by the standard error of the corresponding observational data, which can be acquired by Eq. 8, where N_i represents the sample size of study i.

We established a covariate model to screen for factors that potentially affect PFS, such as age, ECOG status, biomarker status, and study design, including blindness and combination drug information. The missing information in each arm was imputed using the median value of the entire study population. Covariates missing proportions of more than 30% were not tested.

Continuous covariates are described in Eq. 9, and categorical covariates are described in Eq. 10. P p o p = P t y p i c a l × e ( C O V − C O V m e d i a n ) × θ C O V (9) P p o p = P t y p i c a l × ( 1 + θ C O V ) (10)

In Eqs 9, 10, P_pop is the population prediction of the corresponding model parameter, P_typical is the typical value of the model parameter, COV is the reported value of the screened covariate, COV_median is the median value of the screened covariate, and θ_cov is the scaling factor.

Forward and backward elimination approaches were used to test the covariates. The covariates were added stepwise to the model if the objective function value (OFV) decreased by > 2.71 (corresponding to p < 0.1). After defining the full model, the significance of each covariate was tested individually by removing them, one at a time, from the full model. Any covariate that failed to decrease OFV by > 3.84 (corresponding to p < 0.05) was removed from the model.

Goodness-of-fit (GOF) plots, visual predictive check (VPC), and bootstrap were used to identify the model’s performance. GOF plots included the following scatterplots: observation (OBS) vs. individual prediction (IPRED), OBS vs. population prediction (PRED), conditional weighted residual errors (CWRES) vs. PRED, and CWRES vs. time.

The VPC was used to check whether the model was able to reproduce the variability and main trends of the observed data. Typically, in this method, 1,000 datasets based on the final model parameters are modeled using Monte Carlo simulations. The observed data are then compared with the 2.5, 50, and 97.5^th percentiles of the simulated data to assess the predictive capacity of the final model.

Bootstrap, which was performed by 1,000 NONMEM repetitions of the final model, was used to assess the robustness of the model. The bootstrap median parameter values and 95% percentile bootstrap confidence intervals (CIs) were compared with the respective values estimated from the original data. The estimated parameters were considered stable and less affected by any individual study if the values were comparable.

A single-arm meta-analysis was performed to summarize the model parameters of each arm using random-effects models according to different situations. The point estimates of the model parameters and their 90% CI in different scenarios were obtained. Next, the parameter values were randomly selected from the parameter distribution of each scenario using the Monte Carlo simulation method. The process was repeated 10000 times at each time point, and the median value and 90% CI of PFS at the different time points and different scenarios were obtained.

The simulation scenarios included the following: 1) PFS of patients after monotherapy or combination therapy; 2) PFS of patients with different biomarkers; 3) PFS of patients after treatment with different PARPis; 4) PFS of patients in combination therapy with different combination drugs; 5) PFS of platinum-sensitive patients (those who had no progression or relapse within 6 months after chemotherapy) and partially platinum-resistant patients (those who progressed or relapsed within 12 months after chemotherapy).

The modeling process was completed using NONMEM7.3 (Level1.0, ICON Development Solutions, New York, United States). First-order conditional estimation was used to estimate the model parameters, and the bootstrap was executed using Perl-Sons-Nonmem (PsN4.9.0). Model simulation, graph drawing, and meta-analysis were performed using the R software (version 3.6.1, The R Foundation of Statistical Computing, Vienna, Austria). The assessment of the literature quality of the RCTs was completed using RevMan (Version 5.4, Nordic Cochrane Center, Copenhagen, Denmark).

A total of 922 articles were retrieved from the PubMed, EMBASE, and Cochrane Library databases, 206 of which were included in the full-text analysis. A total of eight articles (Kaye et al., 2012); (Ledermann et al., 2014); (Oza et al., 2015); (Mirza et al., 2016); (Coleman et al., 2017); (Pujade-Lauraine et al., 2017); (Mirza et al., 2019); (Liu et al., 2019) met the inclusion criteria (Supplementary Figure S1; Supplementary Material).

Of the eight included studies, there were 1,169 patients with ovarian cancer who received PARPis. The age range of the patients was 53.5-67.0 years (median, 58.1 years). There were 501 patients receiving olaparib treatment, 302 receiving niraparib treatment, and 366 receiving rucaparib. Among them, 669 patients tested positive for BRCA mutations (BRCAm) and a 100 patients tested positive for wild-type BRCA (BRCAwt). 400 patients were HRD-positive. A total of 1,077 patients received monotherapy, and 92 patients received combination therapy. The combined drugs included anti-angiogenic drugs: 21 patients and 43 patients for bevacizumab and cediranib, respectively, and 28 patients received a PARPi combined with a chemotherapy drug. There were 64 partially platinum-resistant ROC patients and 1,105 platinum-sensitive ROC patients. Demographic characteristics are shown in Table 1. Detailed information of each study is shown in Supplementary Table S1 (Supplementary Material).

The overall quality of the included studies was high, and there was no obvious risk of bias. However, it should be noted that the four included studies were open-label trial designs (Supplementary Figure S2; Supplementary Material).

The log-normal hazard function model was finally selected as the basic hazard function model as per the OFV minimization and precision of model parameter estimation.

When covariates were examined, we found that the type of therapy (monotherapy or combination therapy) had a significant impact on the model parameter μ . The final model is expressed in Eq. 11, and the model parameters are listed in Table 2. μ , p o p = 2.96 × ( 1 + θ C o m b o ) where: θ C o m b o = { − 0.215   m o n o t h e r a p y 0   c o m b i n a t i o n   t h e r a p y 11

The GOF plots of the final model show that OBS and PRED, as well as OBS and IPRED, are evenly distributed on both sides of the diagonal, and that the regression trend line and standard line coincide. The CWRES of most points are evenly distributed around the 0 line within 6, and the fitting line of CWRES vs. time and PRED almost coincide with the 0 line, indicating that the model fits the measured data at different time points or at different observations, without obvious bias (Supplementary Figure S3, Supplementary material).

The success rate of the bootstrap resampling of the model was 96.3%. In addition, the model parameter distribution obtained by the bootstrap was close to the estimated model parameter value from the original dataset (Table 2), indicating that the model parameter estimation was relatively robust and less affected by individual studies. The VPC results showed that the 95% CI of model prediction covered most of the measured values of drug efficacy, indicating that the model had good predictive performance (Supplementary Figure S4; Supplementary Material).

Sensitivity analysis showed that after deleting all the data from open-label studies, the estimation of the model parameters was consistent with the original data set, suggesting that the model was more robust and less affected by open-label studies(Table 2).

Based on the final model, we estimated the typical PFS of different types of treatment (monotherapy or combination therapy) (Figure 1). Simulation results showed that the typical PFS of monotherapy and combination therapy was 8.5 months (90% CI: 7.5-10.0 months) and 16 months (90% CI: 14.0-19.5 months), respectively.

In this study, 501 subjects were treated with olaparib, 302 with niraparib, and 366 with rucaparib. The median PFS values of olaparib, niraparib, and rucaparib monotherapy were approximately 9.5, 10.5, and 12.0 months, respectively. For combination therapy, the median PFS values of olaparib, niraparib, and rucaparib were approximately 19.0, 20.0, and 25.0 months, respectively.

A total of 669 patients with BRCA mutations, 100 subjects with wild-type BRCA, and 400 HRD-positive subjects were included in the study. The median PFS of patients with BRCA mutation, wild-type BRCA subjects, and HRD-positive subjects were approximately 11.0, 7.5, and 9.0 months in monotherapy and 23.0, 14.0 and 17.5 months in combination therapy, respectively.

A total of 64 patients had partially platinum-resistant ROC and 1,105 patients had platinum-sensitive ROC. The median PFS of patients with platinum-sensitive ROC was approximately 11.0 and 22.0 months with monotherapy and combination therapy, respectively, and that of patients with partially platinum-resistant ROC was approximately 5.5 and 9.0 months with monotherapy and combination therapy, respectively.

In the combination study, patients treated with combination drugs, including antiangiogenic drugs such as bevacizumab and cediranib, were 28 and 44, respectively, and 20 were treated with chemotherapeutic drugs. The median PFS values in combination with cediranib, bevacizumab, and chemotherapy were approximately 17.0, 12.5 and 19.5 months, respectively. All the data are listed in Figure 2 and Table 3.