Population pharmacokinetics of bevacizumab in cancer patients with external validation

Background Bevacizumab is approved for various cancers. This analysis aimed to comprehensively evaluate bevacizumab pharmacokinetics and the influence of patient variables on bevacizumab pharmacokinetics. Methods Rich and sparse bevacizumab serum concentrations were collected from Phase I through IV studies in early and metastatic cancers. Bevacizumab was given intravenously as single agent or in combination with chemotherapy for single- and multiple-dose schedules. Results Model-building used 8943 bevacizumab concentrations from 1792 patients with colon/colorectal, non-small cell lung, kidney, pancreatic, breast, prostate and brain cancer. Bevacizumab doses ranged from 1 to 20 mg/kg given once every 1, 2 or 3 weeks. A two-compartment model best described the data. The population estimates of clearance (CL), central volume of distribution (V1) and half-life for a typical 70-kg patient were 9.01 mL/h, 2.88 L and 19.6 days. CL and V1 increased with body weight and were higher in males than females by 14 and 18 %, respectively. CL decreased with increasing albumin and decreasing alkaline phosphatase. The final model was externally validated using 1670 concentrations from 146 Japanese patients that were not used for model-building. Mean prediction errors were −2.1, 3.1 and 1.0 % for concentrations, CL and V1, respectively, confirming adequate predictive performance. Conclusions A robust bevacizumab pharmacokinetic model was developed and externally validated, which may be used to simulate bevacizumab exposure to optimize dosing strategies. Asian and non-Asian patients exhibited similar bevacizumab pharmacokinetics. Given the similarity in pharmacokinetics among monoclonal antibodies, this may inform pharmacokinetic studies in different ethnic groups for other therapeutic antibodies. Electronic supplementary material The online version of this article (doi:10.1007/s00280-016-3079-6) contains supplementary material, which is available to authorized users.


Introduction
Bevacizumab (Avastin ® , Genentech Inc.) is a humanized monoclonal immunoglobulin G (IgG) 1 antibody that specifically binds and neutralizes the biological activity of vascular endothelial growth factor A (VEGF-A), a key isoform of VEGF involved in angiogenesis, and a well-characterized pro-angiogenic factor [1]. Bevacizumab causes inhibition of tumor angiogenesis by blocking VEGF-A from binding to its receptors and leads to tumor growth inhibition. Bevacizumab in combination with standard therapy has received marketing authorization for use in the treatment of various cancers including metastatic colorectal cancer (CRC) [2,3], non-small cell lung cancer (NSCLC) 1 3 [4], breast cancer [5], renal cell carcinoma [6], cervical cancer [7] and ovarian cancer [8].
A population pharmacokinetic (PK) model has been previously developed [9]. Bevacizumab PK showed dose linearity within the dose range of 1-20 mg/kg, a slow clearance, a volume of distribution consistent with limited extravascular distribution and a terminal half-life of approximately 20 days. Clearance (CL) and central volume of distribution (V1) increased with body weight and were higher in male patients. CL decreased with increasing albumin and decreasing alkaline phosphatase. There has been no evidence for anti-therapeutic antibodies (ATAs) for bevacizumab in metastatic solid tumors based on the large number of historical clinical studies, and ATA was detected in only 0.6 % of the patients with colon cancer (adjuvant setting) [10].
However, the previous analysis had several limitations. Several important covariates were not evaluated in previous analysis including ethnicity (e.g., Asian vs. non-Asian), indications and baseline VEGF-A. First, bevacizumab has been widely used across ethnic groups (e.g., Asian vs. non-Asian), and supplementary Biologics License Applications have been submitted to health authorities for approval of using bevacizumab for new indications or new combinations based on data from limited ethnic groups while the target population contains much broader ethnic groups. Therefore, it is important to evaluate ethnicity (Asian vs. non-Asian) as a covariate. Second, it has been shown that bevacizumab clearance is 50 % higher and exposure is 50 % lower in gastric cancer as compared to other types of solid tumors [11], making it important to evaluate indication as a covariate. Finally, several studies have shown the predictive value of baseline VEGF-A for bevacizumab treatment effect on progression-free survival and/or overall survival, meaning that only patients with high VEGF-A levels may benefit from bevacizumab treatment, for example in gastric cancer [12] and metastatic breast cancer [13]. Therefore, it is important to evaluate baseline VEGF-A as a covariate. The reason why these important covariates were not evaluated in the previous analysis is likely that these evidences showing the importance of these covariates all appeared after the previous analysis, and therefore, the significance of these covariates may not have been fully realized at that time, and/or the data were unavailable at that time.
Other limitations of the previous analysis include utilization of FO (first-order) instead of FOCE (first-order conditional estimation) algorithm in NONMEM [14], limited number of studies (n = 6), patients (n = 491) and indications (mainly CRC, NSCLC and breast cancer), etc. Therefore, an updated analysis is warranted.
The objectives of the current analysis were to develop a robust population PK model in adult patients with solid tumors and to evaluate the influence of patient variables on bevacizumab PK, which can be used to simulate bevacizumab exposure to optimize bevacizumab dosing strategies.

Patients
Studies in adult cancer patients included in this analysis are summarized in Table 1. Patients with at least one PK sample were evaluated. Serum bevacizumab concentrations were determined at Genentech, Inc., using an enzymelinked immunosorbent assay that used recombinant human VEGF for capture and a goat antibody to human IgG conjugated to horseradish peroxidase for detection. The lowest limit of quantification (LLOQ) was 78 ng/mL in serum [9]. Concentrations below the LLOQ were omitted. The clinically relevant covariates tested included those related to demographics, biochemical tests, concomitant medications and pathophysiological factors ( Table 2). Values of covariates that follow lognormal distribution were  [15] and R 3.0.3 [16]. Several models with various residual error structures and OMEGA matrices were tested to select the optimal base model. The base model included a power function of body size [e.g., total body weight (BWT)] on all PK parameters: where BWT i = baseline BWT of patient i; P i = typical PK parameters of patients with BWT i ; P TV = typical value of PK parameters for patients with BWT of 70 kg; θ P = exponent for the PK parameter P. The base model was evaluated using either theoretical (0.75 for clearances and 1 for volumes of distribution) [17] or fitted values of exponents θ P . The quality of fit was evaluated using a standard model discrimination process including statistical criteria [i.e., minimum of objective function value (OFV)], adequate estimation of the parameters (e.g., relative standard error <50 %) and graphical representations of goodness of fit. The final model was established in a stepwise manner by forward addition followed by backward elimination of parameter-covariate relationships with a significance level of p < 0.01 and p < 0.001, respectively (OFV decrease of 6.63 and 10.83 for one degree of freedom based on Chi-squared distribution, respectively).
The effect of n covariates at baseline on PK parameters was coded using a multiplicative model: where θ i is the typical value of the parameter for patients with a set of covariates i, θ TV is the typical value of the PK parameter for patients having the covariate values equal to the median of the covariate for all patients, and Effect 1,i through Effect n,i are multiplicative factors of the effects for covariate 1 through n, for the set of covariates i. The covariate models for both continuous and categorical covariates were chosen to avoid prediction of negative parameter values.
The multiplicative factor was defined using the power function for continuous covariates: and defined as follows for categorical covariates: if this categorical covariate is equal to 0, then Effect i = 1 if this categorical covariate is not equal to 0, then where Effect i is the multiplicative factor of the covariate effect for covariate i at baseline, Cov i is the covariate value, Cov reference value is the median of the covariate for all patients, and θ eff is the exponent of the power function to be estimated. e θeff was used for categorical covariates to force a positive value.
For patients with missing value for a continuous covariate Cov i , the multiplicative factor of Cov i was calculated as [18]: where θ cov missing is the value of the covariate that was estimated by fitting to the data from all patients with missing information. During the forward addition step, θ cov missing was estimated by fitting to the data for all covariates with missing values. During the backward elimination step, θ cov missing was fixed to the estimates obtained from the forward addition step to stabilize the model. All patients with missing value had the same estimate of θ cov missing .

Model evaluation and sensitivity analysis
The population PK models were evaluated using diagnostic plots [19,20], visual predictive check (VPC) [20, 21], prediction-corrected VPC (pcVPC) [22], bootstrapping [23] and shrinkage [24] assessments. VPC, pcVPC and bootstrapping were all performed using 2000 replicates based on the final model. The relative impact of each covariate included in the final model alone on PK parameters and exposure was explored. Exposure including the trough (C min ) and peak (C max ) concentration was computed at steady state given bevacizumab of 10 mg/kg once every 2 weeks using the final model. The computation was performed using the extreme covariate values (5th and 95th percentiles) and the equations of the final model.

External validation
After the final model was built, data from three Japanese studies (JO18157, JO19901 and JO19907, Tables 1 and 2) became available and were subsequently used for external validation. Predicted bevacizumab concentrations (C IPRED ) for the validation population were obtained using post hoc Bayesian forecasting by fixing the parameters in the structural and variance models to the final estimates. Prediction errors (PEs) as a measure of bias were calculated for each concentration as PE = (C IPRED − C OBS )/C OBS , where C OBS denotes observed concentrations. Root mean squared prediction error (RMSE) as a measure of precision was calculated where n denotes the number of observations. pcVPC approach was used to compare the 95 % prediction interval and OBS.
Predicted PK parameters (P IPRED ) for each patient were calculated based on individual covariate values using the equations in the final model without considering observed concentrations. Post hoc estimates of PK parameters (P EST ) were obtained based on observed concentrations and the final model. PE was calculated as (P IPRED − P EST )/P EST . RMSE was calculated as 1 n (P IPRED − P EST ) 2 , where n denotes the number of patients.

Patients
A total of 8943 bevacizumab serum concentrations from 1792 adult cancer patients in 15 studies were included in the model-building dataset, and 1670 concentrations from 146 adult patients in three Japanese studies were included in the external validation dataset. Studies and patient characteristics are summarized in Tables 1 and 2. Less than 5 % of the samples were below LLOQ, and all of them were pre-dose samples.

Population pharmacokinetic modeling
The optimal base model was a linear two-compartment model with theoretical exponents estimated for clearance (CL), inter-compartment clearance (Q), central (V1) and peripheral (V2) volumes of distribution, full block interindividual variability (IIV) on CL, V1 and V2 with combined additive and proportional residual error. Parameter estimates of the base model are presented in Supplementary Table 1. In the base model, the estimates of typical bevacizumab CL, V1, Q, V2 and terminal half-life values for a 70-kg patient were 9.01 mL/h, 2880 mL, 18.7 mL/h, 2571 mL and 19.6 days. Thirty-eight covariate relationships were evaluated in the forward addition step. After adjusting for total body weight (BWT), CL and V1 were still higher in male patients. CL decreased with increasing albumin (ALBU) and decreasing baseline alkaline phosphatase (BALP). CL was also lower in patients treated with interferon alpha (Supplementary Fig. 1). No covariate was excluded during the backward elimination step (p < 0.001).
Parameter estimates of the final model are summarized in Table 3. Bevacizumab CL and V1 for the patient i were described as follows (ALBU = 41.8 and BALP = 76.3 if missing):

Model evaluation and sensitivity analysis
Goodness-of-fit plots showed good agreement between predicted and observed bevacizumab concentrations with no apparent bias in residual ( Supplementary Fig. 2). The pcVPC result for the final model is presented in Fig. 1. Overall, the 2.5th, 50th and 97.5th percentiles of observed concentrations were within the predicted 95 % confidence interval (CI) of these percentiles, suggesting accurate model fitting across a wide range of dosing regimens and time courses. Bootstrapping resulted in median parameter estimates and 95 % CIs similar to the estimates from the original dataset, indicating that the final model provided good precision for parameter estimation.
The impact of the variation for a single covariate included in the final model on steady-state exposure   (Fig. 2a, b), CL (Fig. 2c) and V1 (Fig. 2d) and steady-state exposure (Fig. 2)

External validation
Over 95 % of prediction-corrected observations fell within the 95 % prediction interval (PI) (Fig. 3). CL and V1 calculated based on the equations in the final model (P IPRED ) were similar to those estimated based on observed concentrations (P EST ) (Fig. 4a, b) Post hoc Bayesian estimates of CL and V1 in the modelbuilding population (mostly non-Asian patients) and external validation population (Japanese patients only) were also similar after normalization by individual covariate values that were included in the final model (Fig. 4c, d). The normalization was done by dividing the post hoc Bayesian estimates of CL and V1 by the individual covariate values in the form that appeared in the equations of the final model.

Discussion
This analysis is a comprehensive PK evaluation of bevacizumab in adult cancer patients in Phases I-IV studies as a single agent or in combination with chemotherapy for both single-and multiple-dose administration with both rich and sparse bevacizumab serum concentration data. A robust population PK model was built based on a large PK population of 1792 patients from 15 studies and then externally validated using data from 146 Japanese patients in three independent studies. This model consolidated all bevacizumab PK data in one model, can timely support simulations and decision making when needed, can help develop consistent pharmacokinetic messages of bevacizumab for investigators and health authorities given that multiple PK models have been developed for bevacizumab and contained inconsistent messages, and can support future studies of bevacizumab in other indications. As mentioned in the "Introduction," many important covariates that were not evaluated in the previous analysis (e.g., Asian vs.  . Red vertical lines represent the "base" defined as the exposure or PK parameter estimate of a typical patient, i.e., a 70-kg female patient with albumin of 39 g/L and baseline alkaline phosphatase of 109 U/L without interferon alpha treatment. The dark blue shaded curve at the bottom with value at each end shows the 5th to 95th percentile range of exposure or PK parameter estimate across the entire population. Each light blue shaded bar represents the influence of a single covariate on the steady-state exposure after repeated bevacizumab dose of 10 mg/kg once every 2 weeks or on the PK parameter. The label at left end of the bar represents the covariate being evaluated. The upper and lower values for each covariate capture 90 % of the plausible range in the population. The length of each bar describes the potential impact of that particular covariate on bevacizumab steady-state exposure or PK parameters, with the percentage value in the parentheses at each end representing the percent change from the "base." The most influential covariate is at the bottom of the plot for each exposure metric or PK parameter 1 3 Typical population PK parameter estimates were similar as previously published [9]. The low IIV of 29 and 18.3 % observed for CL and V1 was typical for antibody drugs [25]. The pcVPC demonstrated adequate fit and predictive performance of the final model (Fig. 1). The median prediction (blue band) may appear to be slightly below the median observation (blue line) beyond day 112, suggesting a possible tendency of under-prediction. However, this tendency is likely irrelevant given (1) the small degree of under-prediction, (2) the sparseness of data beyond day 112, (3) the good predictive performance for 2.5th and 97.5th percentiles (red) across all time points, as well as (4) the complexity and heterogeneity of the data. The external validation (Figs. 3, 4) demonstrated good predictive performance of the final model with no apparent systemic bias and the similarity in bevacizumab PK between Asian and non-Asian adult cancer patients. Although there may appear to be a tendency of over-prediction of the variability (Fig. 3), this tendency is likely irrelevant because the model was built based on more heterogeneous data (15 studies over a decade across various ethnic groups) while the validation data were more homogeneous (three studies over a few years in Asian patients only).
Factors significantly associated with bevacizumab PK were similar as previously published [9]: CL and V1 increased with BWT and were higher in males, and CL decreased with increasing albumin and decreasing BALP.
It is well known that CL of other IgG antibodies is faster in patients with lower serum albumin levels [25], likely due to two reasons. First, the level of albumin correlates with disease status. Second, the recycling of albumin and IgG is both mediated by FcRn (neonatal Fc receptor) [25], and therefore, albumin levels may reflect the abundance and efficiency of FcRn. The effect of BALP on bevacizumab CL is likely because BALP is an indicator of disease burden, such as liver or bone metastases. CL was found to be 15.6 % lower in patients treated with interferon alpha. However, this effect was within the overall PK variability and therefore may be clinically irrelevant.
Similar to the previous analysis [9], tumor burden was not included in the final model in this analysis. Among solid tumors, tumor burden is usually an indicator of disease severity and health status. It is usually defined as the sum of longest diameters of target lesions under RECIST (Response Evaluation Criteria in Solid Tumors) criteria for systematic tumors and under other criteria for other tumors (e.g., brain tumors). Inclusion of tumor burden in bevacizumab PK model may not be crucial. First, tumor burden as an indicator of disease burden and health status could already be represented by albumin and BALP in the model. Second, tumor burden as a source of VEGF-A (target of bevacizumab) is irrelevant for bevacizumab PK because bevacizumab molar concentration is thousands of times higher than that of VEGF-A [10], and there has been no evidence of target-mediated drug disposition (TMDD) for bevacizumab [10]. Third, in previous analyses, tumor  Fig. 4c, d, data points with CL < 3 mL/h (n = 2) or V1 < 1500 mL (n = 1) are not displayed 1 3 burden alone showed relatively low impact on bevacizumab exposure in the sensitivity analysis (similar to Fig. 2, data not published). Finally, the final model demonstrated adequate fitting and superior predictive performance without incorporating tumor burden.
On the other hand, three factors made it impossible to test baseline tumor burden as a covariate in this analysis. First, tumor response criteria were inconsistent across these 15 studies that were conducted across a time span of over a decade. Several different versions of RECIST and other criteria (e.g., Macdonald criteria for glioblastoma in BO21990) were used. Second, the methods used to measure tumor burden were inconsistent across studies, such as CT (computerized tomography) scans and MRI (magnetic resonance imaging). Finally, unit of length (mm) and area (mm 2 ) both exist in tumor burden data, which cannot be converted to each other. In fact, inclusion of tumor burden in the model would greatly reduce the applicability of the model due to the continuous advancement in tumor response criteria and measurement methods, and due to different tumor response criteria and measurement methods across cancer types, for example RECIST version 1.0 versus version 1.1, RECIST criteria versus Macdonald criteria or RANO (Response Assessment in Neuro-Oncology) criteria, CT scans versus MRI.
In conclusion, a robust population PK model for bevacizumab in adult cancer patients was built and externally validated, which may be used to simulate concentrationtime profile in adult cancer patients in future studies [11,26]. Baseline body weight, albumin, alkaline phosphatase and gender were the covariates with the greatest influence on bevacizumab CL and V1, supporting body weight-based dosing of bevacizumab. No difference in bevacizumab PK was observed between Asian and non-Asian patients. Given the similarity in PK among many monoclonal antibodies, this may inform PK studies in different ethnic groups (e.g., Asian vs. non-Asian) for other therapeutic antibodies without TMDD and significant race-dependent target expression.
Author contributions All authors have contributed substantially to conception and design of the analysis, drafting or revising the paper as well as giving final approval for submission.

Compliance with ethical standards
Conflict of interest Kelong Han, Angelica Quartino, Sandhya Girish, David E. Allison and Jin Jin receive salary from Genentech and hold stocks of Roche Pharmaceuticals; in addition, Jin Jin also holds stock in Eli Lilly; Thomas Peyret, Mathilde Marchand and Nathalie H. Gosselin receive salary from Pharsight Consulting Services; all authors declare: no financial relationships with any organizations that might have an interest in the submitted work in the previous 3 years; no other relationships or activities that could appear to have influenced the submitted work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.