A pharmacokinetic binding model for bevacizumab and VEGF165 in colorectal cancer patients

Purpose To characterize the population pharmacokinetics of bevacizumab, its binding properties to VEGF165 and the effect of demographic data and VEGF-A polymorphisms on the interplay between bevacizumab serum pharmacokinetics and VEGF165 serum concentrations in patients with colorectal cancer stage IV. Methods Bevacizumab and VEGF165 data were collected from 19 adult patients with metastatic colorectal cancer enrolled in an observational clinical study. Bevacizumab was administered with one of the following combinations: 5-FU/Leucovorin/Irinotecan, 5-FU/Leucovorin/Oxaliplatin, Capecitabine/Irinotecan at doses ranging from 5 to 10 mg/kg every 2 or 3 weeks. Data analysis was performed using nonlinear mixed-effects modeling implemented in NONMEM 7.3. Results A target-mediated drug disposition model adequately described bevacizumab concentration changes over time and its binding characteristics to VEGF165. The estimated clearance of bevacizumab was 0.18 L/day, the free VEGF165 levels at baseline were 212 ng/L, and the elimination rate constant of free VEGF165 was 0.401 day−1. Body weight was allometrically included in all PK parameters. Conclusion The final model adequately described the pre- and post-dose concentrations of total bevacizumab and free VEGF165 in patients with colorectal cancer. Model parameters were consistent with those previously reported for patients with solid tumors. Correlations between the binding affinity of bevacizumab and the VEGF-2578C/A and VEGF-634G/C polymorphisms were noticed. Electronic supplementary material The online version of this article (doi:10.1007/s00280-015-2701-3) contains supplementary material, which is available to authorized users.

There are five major isoforms of the human VEGF-A glycoprotein (VEGF 121 , VEGF 145 , VEGF 165 , VEGF 189 , VEGF 206 ) produced by the alternative splicing of the VEGF-A gene. All isoforms differ in the amino acid length and the binding ability to heparin and heparin sulfate proteoglycans found on the cell surface or in the extracellular matrix. VEGF 145 , VEGF 189 and VEGF 206 are tightly bound to the cell surface and the extracellular matrix, VEGF 121 is freely diffusible, and VEGF 165 is moderately diffusible. VEGF 165 is the most abundant and potent stimulator of angiogenesis compared to the other isoforms. It is overexpressed in several cancer types and is related to the ability of a tumor to grow, invade and spread [7,8].
Several clinical trials in CRC patients indicate an inverse relationship between the concentration of VEGF-A in serum or tumor tissue and the clinical outcome [9][10][11]. Furthermore, certain VEGF-A polymorphisms have been associated with altered VEGF-A production or promoter activity, causing variability in response to treatment [12][13][14]. The VEGF-2578CC, VEGF-1154GG and VEGF-634CC genotypes have been correlated with higher VEGF-A production than other genotypes for VEGF-2578C/A, VEGF-1154G/A and VEGF-634G/C polymorphisms [13,15]. Thus, free serum VEGF-A levels could serve as a useful biomarker for reflecting the anti-angiogenic activity of bevacizumab and for indicating when dose adjustments are needed to achieve sufficient VEGF-A blockade.
Monoclonal antibodies exhibit more complex concentration-time [pharmacokinetic (PK)] and effect-concentration [pharmacodynamic (PD)] characteristics than small molecules. The main elimination pathway is proteolytic catabolism throughout the body (linear, nonspecific clearance) and not hepatic metabolism or renal excretion. However, the disposition of monoclonal antibodies is often influenced by their high affinity binding to their molecular targets and subsequent degradation of the monoclonal antibody-target complexes via fluid phase or receptormediated endocytosis. This phenomenon is known as target-mediated drug disposition (TMDD), and it usually leads to nonlinear, saturable distribution and elimination [16,17]. In recent years, population PKPD modeling has proven to be a useful tool to elucidate the dose-response relationships, to explain inter-individual variability in the observed drug exposure or response and to guide dose selection to achieve the optimal benefit-risk ratio for a given anticancer treatment [18][19][20]. TMDD modeling has been generally used to describe the dynamics of the drug-target interaction, as it can provide valuable information on the mechanism of the underlying PKPD relationship [17,21,22]. The full TMDD model is a complex system of differential equations describing the concentrations of the free drug, the free target and the drug-target complex. Model parameters may not be identifiable when limited experimental data on target, total and free drug concentrations are available. In that case, simplifications of the full TMDD model with fewer parameters can be applied such as the quasi-equilibrium, quasi-steady-state (QSS) and Michaelis-Menten approximation [23,24]. The full TMDD model is based on the following assumptions: (1) the drug-target binding is a simple, one-to-one binding process, (2) the drug can only bind to its specific molecular target, (3) the drug-target binding takes place in the central and not in the peripheral or depot compartments, (4) only the free drug can distribute into the peripheral compartment, (5) the drug-target complex is totally eliminated, and (6) the free target production and degradation rates are constant and independent of the drug or target concentrations. The main assumption of the QSS approximation, which differentiates it from the full TMDD model or other TMDD approximations, is that the change in the drug-target complex concentration is negligible on the time scale of the other system processes. The free drug, the free target and the drug-target complex are assumed to be in a QSS. The QSS approximation is usually preferred when only total drug and target concentrations are available and is applied for drugs with fast binding, dissociation and internalization rates [22][23][24].
The aim of the present study was to describe the PK of bevacizumab, its binding properties to VEGF 165 and the effect of individual patient characteristics on the relationship between bevacizumab and VEGF 165 in patients with stage IV CRC using a population modeling approach. To the best of our knowledge, a limited number of population PK studies for bevacizumab have been conducted so far [25,26]. Moreover, little quantitative pharmacology research has been performed to establish the impact of VEGF-A polymorphisms and other relevant covariates on treatment response [27][28][29]. Therefore, it is anticipated that the development of a binding model for bevacizumab in patients with metastatic CRC would yield a better mechanistic understanding of the biological system, as it could characterize the in vivo interaction of the drug with its soluble target, VEGF 165 . It would further elucidate the PK characteristics of the drug, its effect on the free target time course and provide more information on target affinity and the influence of covariates (e.g., genetics).

Patients and study design
Nineteen subjects were enrolled in this observational (non interventional) study. Eligible patients had CRC stage IV documented by histology (biopsy through colonoscopy or surgery), CT scans (thoracic, abdominal and brain scans), bone scintiscans and the measurement of cancer biomarkers (CEA and CA 19.9); Eastern Cooperative Oncology Group (ECOG) performance status ≤2; no chemotherapy or radiotherapy within 3 months of first study treatment; were at least 18 years of age. Prior to enrollment, all patients had undergone a complete history and physical examination. Standard hematologic and biochemical laboratory tests had also been conducted to assess an adequate bone marrow, liver and renal function as defined by: white blood cells count (WBC) ≥ 2500/μL, absolute neutrophil count (ANC) ≥ 1500/μL, platelets ≥ 75,000/μL, serum creatinine ≤ 1.5 × upper limit of normal (ULN), urine protein ≤ 2 g/24 h, total serum bilirubin ≤ 3 mg/dL, AST (SGOT)/ALT (SGPT) ≤ 2 × ULN. The research was conducted in accordance with the Declaration of Helsinki, and all the appropriate approvals were obtained by the relevant Ethics Committees. Signed informed consent was provided by all participants before the initiation of the study.
The current study was carried out in three Greek oncology clinics, at the Department of Medical Oncology of Agii Anargiri Cancer Hospital, at the Department of Medical Oncology of Metropolitan Hospital in Athens and at the Division of Oncology of Patras University Hospital in Rio.
Bevacizumab (Avastin ® , Roche Registration Ltd.) was administered as an intravenous infusion at a dose of 5 mg/ kg, in combination with 5-Fluorouracil/Leucovorin/Irinotecan (BEV-FOLFIRI treatment) or 5-Fluorouracil/Leucovorin/Oxaliplatin (BEV-FOLFOX treatment), in 2-week cycles or at a dose of 7.5 mg/kg together with Capecitabine/Irinotecan (BEV-CAPIRI treatment) in 3-week cycles. One patient received the BEV-FOLFIRI treatment at a dose of 10 mg/kg in 2-week cycles. The duration of first infusion was 90 min. If no infusion-related symptoms were observed, subsequent infusions were given over 60 or 30 min.
Subjects on BEV-FOLFIRI or BEV-FOLFOX treatment received initially six cycles of bevacizumab and those who responded to treatment, as defined by Response Evaluation Criteria in Solid Tumors [30], continued with the same treatment for other six cycles. The responders, after the completion of 12 cycles in total, continued with bevacizumab monotherapy till disease progression. In that case, they were administered bevacizumab at a dose of 5 mg/kg every 2 weeks. Subjects on BEV-CAPIRI treatment received initially three or four cycles of bevacizumab (the number of cycles depends on the lines of treatment the patients had undergone before the initiation of the current treatment). After the completion of eight or nine cycles, the responders continued either with the same treatment or bevacizumab monotherapy upon development of disease progression at the recommended dose (7.5 mg/kg every 3 weeks).
Pre-and post-dose (after the end of infusion) concentrations of total bevacizumab (free and bound to one molecule of VEGF 165 ) and free VEGF 165 (unbound to bevacizumab) were measured in serum during several cycles of treatment. Two blood samples (one pre-dose and one post-dose) were drawn on day 1 of the following cycles: (1) 3, 6, 8, 12, 18 and 24 (patients on BEV-FOLFIRI or BEV-FOLFOX treatment), (2) 2, 4, 5, 8 and 11 (patients on BEV-CAPIRI treatment). One pre-dose blood sample was also collected on day 1 of cycle 1, intended only for free VEGF 165 analysis. The sample collection schedule for total bevacizumab and free VEGF 165 measurements is shown in Supplemental Fig. 1. The genotypes for VEGF-2578C/A, VEGF-1154G/A and VEGF-634G/C single-nucleotide polymorphisms (SNPs) were also determined in blood.

Measurement of total bevacizumab in serum
Blood samples were collected in serum separator tubes and were allowed to clot for 30 min. After centrifugation at 1000×g for 20 min, the serum was removed and stored in aliquots at ≤−20 °C until analysis.
The concentration of total (free and bound to one molecule of VEGF 165 ) bevacizumab in serum was measured using a previously published enzyme-linked immunosorbent assay (ELISA), where the detection limit was 0.033 mg/L and the range of linearity was between 5 and 75 mg/L with precision 5.6 % [expressed as coefficient of variation (CV) percentage]. Standards of 0.24, 0.47, 0.94, 1.88, 3.75, 7.5, 15 and 30 mg/L were used to generate the standard curve, which are well above the detection limit of the assay and within the range of linearity [31].
Microtiter Nunc Maxisorp 96-well plates were coated with recombinant human VEGF 165 (R&D Systems ® Europe) at a concentration of 0.15 mg/L in carbonate-bicarbonate buffer (1 M, pH 9.6) overnight at 4 °C (100 μL/well). After washing four times with phosphate-buffered saline (PBS) containing 0.05 % Tween 20, the wells were blocked with PBS containing 1 % BSA (200 μL/well) and were incubated for 2 h at room temperature. Afterward, the plates were washed and 100 μL of 1:100 diluted standards and samples in 1 % PBS-BSA was added and were incubated for 1 h at 37 °C in an incubator shaker. Then, the plates were washed again, and 100 μL of peroxidase-conjugated goat antihuman IgG specific for Fc fragment (AbD Serotec ® , A Bio-Rad Company) diluted in 1 % PBS-BSA was added to each well. After 1-h incubation at room temperature followed by washing, 100 μL OPD (Sigma-Aldrich) was added and the reaction was allowed to develop at room temperature in the dark. The color reaction was stopped with the addition of sulfuric acid (2 M, 50 μL/well).
The optical density was measured at 450 nm with a correction at 650 nm using an ELISA plate reader (Ther-moMax, Molecular Devices). Duplicate readings for 1:100 diluted standards and samples were performed.
The best fit line of the standard curve was determined by regression analysis using OriginPro 8.0 software (Origin-Lab ® Corporation). The concentrations read from the standard curve were multiplied by the dilution factor.

Measurement of free VEGF 165 in serum
Blood samples were collected in serum separator tubes and were allowed to clot for 30 min. After centrifugation at 1000×g for 20 min, the serum was removed and stored in aliquots at ≤−20 °C until analysis.
The concentration of free VEGF 165 (unbound to bevacizumab) in serum was measured by a commercially available ELISA kit for VEGF 165 (Quantikine ® human VEGF, R&D Systems ® Europe). The detection limit of the assay was 9 ng/L, and the precision was 6.7 % (CV %) [32]. According to the manufacturer, this ELISA assay has not been tested yet for interference with the detection of free or total (free and bound to bevacizumab) VEGF 165 in the presence of bevacizumab. To confirm the hypothesis that it can only discriminate and quantitate free VEGF 165 , we measured VEGF 165 concentrations in samples after the addition of increasing concentrations of bevacizumab. VEGF 165 standards (1000 and 250 ng/L, respectively) were mixed with increasing VEGF 165 -to-bevacizumab molar ratios of 1:0, 1:0.1, 1:1 and 1:1000.
The assay procedure is briefly described below. Plates pre-coated with a mouse anti-VEGF antibody were used to capture VEGF 165 in standards or samples. Any unbound proteins were washed off and a peroxidase-conjugated polyclonal antibody specific for VEGF 165 was added. Then, the plates were washed again and tetramethylbenzidine substrate solution was added. A blue color was developed in proportion to the amount of VEGF 165 present in the ELISA samples. Color development was stopped with the addition of sulfuric acid.
The optical density was measured at 450 nm with a correction at 550 nm using an ELISA plate reader (ELx800™, BioTek Instruments). All standards' and samples' readings were performed in duplicate.
A standard curve was generated with VEGF 165 concentrations ranging from 31.2 to 2000 ng/L. The best fit line was determined by regression analysis using OriginPro 8.0 software (OriginLab ® Corporation).

VEGF genotyping
Genomic DNA was isolated from blood (3 mL) using the Gentra Puregene Blood kit (QIAGEN). DNA concentrations were determined by measuring the optical density at 260 nm with a UV-Vis spectrophotometer (NanoDrop 2000, Thermo Fisher Scientific). DNA purity, which is indicated by the ratio of optical density at 260 and 280 nm, was 1.7-1.9.

Model development
The nonlinear mixed-effects modeling software NON-MEM ® 7.3 (Icon Development Solutions, Ellicott City, MD, USA) was used in the data analysis. All population parameter estimates were obtained with the first-order conditional estimation method with interaction. The graphical representation of the data and model diagnostics was performed with the software tool Xpose (version 4.3.5) [35]. The PsN toolkit (version 4.2.0) was used for the implementation of computer-intensive statistical methods (stepwise covariate model building and randomization test) [36]. Typical concentration-time profiles for total bevacizumab, total and free VEGF 165 were generated using Berkeley Madonna (version 8.3.18, Kagi Shareware, Berkeley, CA, USA).
Inter-individual variability in model parameters was tested assuming a log-normal distribution described by an exponential model (Eq. 1), where P i represents the parameter estimate for the ith individual, P p the typical parameter estimate in the population and η i,p the random variable for the ith individual from a normal distribution with a mean of zero and an estimated variance of ω 2 . Proportional, additive and combined error models on normal scale or additive and combined error models on log-transformed scale were explored to describe the unexplained residual variability. The magnitude of inter-individual and residual variability was expressed as CV %.
The likelihood ratio test was used to assess whether the difference in the objective function (ΔOFV) between different (sub)models (assumed to be χ 2 distributed) was statistically significant. The evaluation of basic goodnessof-fit plots, shrinkage [37] and parameter uncertainty were also taken into account for model discrimination.
Plots of individual empirical Bayes (post hoc) estimates of the PK and VEGF-related parameters versus covariates were explored for the identification of potential parameter-covariate relationships. The choice of the covariate model was based on a stepwise covariate model procedure, a PsN feature that implements forward selection and backward elimination of covariates to a model according to statistical criteria [36]. A decrease in OFV of 3.84 points for forward inclusion (p < 0.05, for one parameter difference) and an increase of 6.64 points for backward deletion (p < 0.01, for one parameter difference) were considered significant.
The following covariates were tested for significance on model parameters: actual body weight, age, gender and VEGF-A SNPs (−2578C/A, −1154G/A and −634G/C). Linear and power parameterizations were considered for the continuous covariates (Eqs. 2, 3, respectively), whereas the categorical covariates were tested in linear equations (Eq. 4). The equations are shown below: In these equations, Cov represents the covariate and Cov med. the median value of the covariate in the study population. θ is the fractional change in the population parameter for an individual with a covariate value different from the median value (Eq. 2) or a reference value (Eq. 4). k stands for the exponential scaling factor and in case of body weight, it was either estimated or fixed to a certain value (0.75 for clearance and 1 for volume parameters) [38]. (1) A randomization test was performed to calculate the probability of identifying falsely significant covariates [39]. A total of 200 new datasets were generated by shuffling 200 times the sequence of the covariate values between individuals in the randomization column. A base model (without a parameter-covariate relationship) was fitted to the original dataset. Then, a full model (with a parameter-covariate relationship) was fitted to the original and the randomized datasets, and the ΔOFVs from the base model were computed. Based on the distribution of the difference in OFV between the two models for each of the 200 datasets, the actual drop in OFV required to reach the 5 % significance level was calculated.
Data analysis was performed in two separate steps. As a first step, only total bevacizumab concentration-time (PK) data were included in the analysis and a population PK model for bevacizumab was developed to gain information on the PK profile of the drug in the current study population. A previously published model by Lu et al. [26], which describes bevacizumab data in a similar study population in terms of demographic and study characteristics, was used as a reference to evaluate whether the present PK data are in reasonable agreement with the previous observations. In the second step of data analysis, a simultaneous fit of total bevacizumab and free VEGF 165 concentration-time data was performed to develop a binding model for bevacizumab using a TMDD modeling approach. The simultaneous data analysis allowed for the exploration of the influence of VEGF 165 on bevacizumab disposition and eliminated the effect of potential shrinkage in PK parameters on the estimation of VEGF-related parameters.

PK model
Total bevacizumab concentrations in serum were expressed in mg/L.
In the first step of data analysis, one-and two-compartment models with first-order elimination were investigated to describe bevacizumab concentration changes over time.
The effect of body weight, age and gender was tested on all PK parameters taking into account prior knowledge on the highly influential covariates [25,26] and the clinical relevance of the available covariates. In addition, the influence of VEGF-A SNPs (−2578C/A, −1154G/A and −634G/C) was only explored on clearance parameters, as these polymorphisms are associated with VEGF production or promoter activity [12][13][14].
Inter-individual variability was investigated in all model parameters, and an error model was developed to describe the unexplained residual variability.
All parameters of the developed PK model were estimated, and the PK parameter estimates reported by Lu et al. [26] were only used as initial estimates. The ability of the proposed model to yield plausible parameter estimates was tested by comparing them with the parameter estimates obtained by the model of Lu et al., which predicts the concentration-time profile of bevacizumab in a larger study population.

TMDD (binding) model
Total bevacizumab and free VEGF 165 concentrations in serum were expressed in nM. Modeling was performed on log-transformed data.
In the second step of data analysis, the structural and covariate form of the developed PK model for bevacizumab was linked with a QSS TMDD model to simultaneously describe total bevacizumab and free VEGF 165 concentration-time profiles.
The QSS equations [22][23][24] for the concentrations of total bevacizumab (C tot ), total VEGF 165 (R tot ), free bevacizumab (C), bevacizumab-VEGF 165 complex (RC) and free VEGF 165 (R) are shown below: In these equations, In(t) is the infusion rate; A 2 is the amount of the free bevacizumab in the peripheral compartment; V 1 and V 2 are the volumes of distribution of the central and peripheral compartment, respectively; CL is the elimination clearance of free bevacizumab; CL RC is the elimination clearance of the bevacizumab-VEGF 165 complex; Q is the intercompartmental clearance. k in and k out are the production (zero-order) and elimination (first-order) rate constants of free VEGF 165 . K ss is the steady-state constant that defines the QSS among free bevacizumab, free VEGF 165 and bevacizumab-VEGF 165 complex and is described by the following equation: k int , k off and k on are the elimination (first-order), dissociation (first-order) and binding (second-order) rate constants of the bevacizumab-VEGF 165 complex. VEGF-A SNPs (−2578C/A, −1154G/A and −634G/C) were only tested for significance on VEGF-related parameters.
Inter-individual variability was investigated in all PK and VEGF-related parameters, and an error model was developed to describe the unexplained residual variability.
All parameters of the developed TMDD model were estimated, and the population PK parameter estimates from the first step of data analysis were only used as initial estimates.

Model evaluation
The predictive performance of the developed models for bevacizumab (PK and TMDD model) was evaluated using prediction-corrected visual predictive checks (pcVPCs) [40,41]. A total of 1000 simulated datasets were generated, and the 95 % confidence intervals for the median, 10th and 90th percentiles of the simulated data were compared to the median, 10th and 90th percentiles of the observed data. The uncertainty in model parameters was expressed by the relative standard error (RSE) obtained by the covariance step in NONMEM.

Patients and samples
The current analysis was based on 86 total bevacizumab and 93 free VEGF 165 serum concentrations collected from 19 adult patients with stage IV CRC. None of the measured concentrations was below the limit of quantification. Blood samples (median 4, 2-10 per patient) were collected before bevacizumab administration and immediately after the end of infusion.
Ten patients received bevacizumab in combination with FOLFIRI or FOLFOX treatment, whereas nine patients were administered bevacizumab together with CAPIRI treatment. Only three patients switched over to bevacizumab monotherapy (one patient from each treatment group). The patient on BEV-CAPIRI treatment had received nine cycles of bevacizumab before changing to bevacizumab monotherapy.
A summary of the patient and study characteristics is shown in Table 1.

VEGF assay
The ability of the Quantikine ® human VEGF ELISA assay [32] to detect free or total VEGF 165 in the presence of bevacizumab was tested in vitro by mixing VEGF 165 standards (1000 and 250 ng/L) with increasing bevacizumab concentrations. The measured VEGF 165 concentrations were equal to zero when bevacizumab was in excess confirming that only free VEGF 165 can be detected by this VEGF ELISA assay, otherwise the measured VEGF 165 concentrations should be identical to the concentrations of VEGF 165 standards (1000 and 250 ng/L). Mean (RSE %) VEGF 165 concentrations for VEGF 165 to bevacizumab 1:0, 1:0.1, 1:1 and 1:1000 molar ratios are shown in Table 2.

PK model
A two-compartment model with first-order elimination was found to describe bevacizumab concentration changes over time better than a one-compartment model (ΔOFV = −18.9). The selection of the final structural model was further supported by basic goodness-of-fit plots and a previously published model for bevacizumab [26]. Goodness-of-fit plots showed no apparent bias in residuals over time or across population-predicted bevacizumab concentrations (Supplemental Fig. 2). V 1 and V 2 were found to be 3.14 and 2.63 L, respectively. Bevacizumab CL was 0.17 L/day, and Q was 0.36 L/ day.
Inter-individual variability could be identified for CL (23 %) and V 1 (15 %). The residual variability was explained by a proportional error model (24 %).
In the stepwise covariate analysis, the effect of several covariates (body weight, age, gender and VEGF-A SNPs) was explored on the PK parameters. None of the available covariates provided an increase of 6.64 points in OFV for backward deletion (p < 0.01, for one parameter difference). Body weight was incorporated with an allometric function in all clearance and volume parameters to ensure model stability considering the strong biological covariate relationship [38,42] and prior information from studies in patients with solid tumors receiving bevacizumab [25,26].
Model parameters could be estimated with acceptable precision (RSE < 25 %) except for V 2 (RSE 50 %), Q (RSE 146 %) and inter-individual variability in V 1 (RSE 38 %). The parameter estimates from the PK model and their RSEs are represented in Table 3.

TMDD (binding) model
A QSS TMDD model was applied to explore the in vivo interaction between bevacizumab and its soluble target, VEGF 165 . In this model, free bevacizumab is eliminated from the central compartment through two pathways: (1) degradation (linear CL) and (2) binding to VEGF 165 and subsequent degradation of the bevacizumab-VEGF 165  complex (linear CL RC ). The structure of the binding model is shown in Fig. 1. The developed model adequately described the time course of total bevacizumab and free VEGF 165 serum concentrations in the current study population. The estimated CL of free bevacizumab was 0.18 L/day, and Q was 1.38 L/ day. The information in the present data was not sufficient to allow for a separate estimation of the elimination clearance of the bevacizumab-VEGF 165 complex, which was therefore set equal to the CL of the free bevacizumab (CL RC = CL). Baseline VEGF 165 (free VEGF 165 at time 0, BM 0 ) was 0.0053 nM corresponding to 212 ng/L (assuming a 1:1 molecular interaction), and k out was 0.401 day −1 . K ss was found to be 267 nM. k in , which is defined as the typical value of BM 0 times the typical value of k out (BM 0,P × k out,P ), was 85 ng/L/day.
Inter-individual variability was assigned to the following parameters: CL, V 1 and BM 0 . The residual variability in total bevacizumab and free VEGF 165 concentrations was explained by an additive error model on log-transformed data.
No direct relationship between VEGF-A SNPs and BM 0 , k out or K ss was identified. However, patients with VEGF-2578AA (homozygous), VEGF-634CC (homozygous) and VEGF-634GC (heterozygous) genotypes seem to have a larger K ss . The effect of VEGF-634G/C polymorphism on K ss becomes significant (ΔOFV = −8.6) only after the inclusion of VEGF-2578C/A polymorphism in K ss (ΔOFV = −6). A plot of ΔOFV versus individuals indicated that one specific individual was driving the relationship between VEGF-2578C/A and K ss (data not shown). A randomization test was then performed to explore whether the inclusion of these covariate relationships was falsely significant. It was concluded that the observed trends between VEGF-2578C/A and K ss or VEGF-634G/C and K ss did not reach statistical significance as the required drop in OFV after the inclusion of VEGF-2578C/A or VEGF-634G/C was 6.8 points or 6.7 points, respectively (p < 0.05, forward inclusion of one parameter).
All the PK and VEGF-related parameter estimates from the binding model and their RSEs are represented in Table 3. The RSEs of the fixed and random effects remained well below 25 % except for the inter-individual variability in V 1 (RSE 36 %), indicating that the parameters could be estimated with acceptable precision. A low ε-shrinkage (10 %) suggests that the model diagnostics are reliable. However, η-shrinkage was relatively high for V 1 (34 %) but <22 % for CL and BM 0 .

Evaluation of the TMDD (binding) model
The TMDD model had a good predictive performance, as expressed by the results of the VPC depicted in Fig. 2. It can adequately describe the time course of total  (14) K ss (nM) 267 (22) Prop. error bev (%) 24 (13) 28 (15) Prop. error VEGF165 (%) 32 (14) Fig. 1 Structure of the binding model for bevacizumab-VEGF 165 interaction. The approximation CL RC = CL was used for purposes of model fitting bevacizumab and free VEGF 165 concentrations as well as the observed variability in the study population.

Simulations
The concentration-time profiles of total bevacizumab, total and free VEGF 165 for a typical patient of 70 kg receiving bevacizumab either at a dose of 5 mg/kg every 2 weeks or 7.5 mg/kg every 3 weeks (the most frequent dosing regimens in the study population) are shown in Fig. 3. A significant drop in the free VEGF 165 levels was observed upon administration of the first dose (73 % for the lower dose and 80 % for the higher dose) followed by a less pronounced decline on subsequent doses attributed to the reversible formation and accumulation of bevacizumab-VEGF 165 complexes (Fig. 3a, b, respectively). The total VEGF 165 concentrations increased over time, in both dosing regimens, up to a level where no more complexes could be formed. The extent of total VEGF 165 accumulation seemed less at the higher dose compared to the lower dose, while the decline in free VEGF 165 levels was slightly increased (54 % compared to 42 % for the lower dose, Fig. 3c).

Discussion
In this study, a binding model for bevacizumab was developed to characterize its PK behavior, to describe its binding properties to VEGF 165 and to assess the influence of relevant covariates (e.g., genetics) on the relationship between bevacizumab and VEGF 165 in adult patients with stage IV CRC. Sparse bevacizumab and VEGF 165 data were collected during routine clinical practice from 19 adult patients following bevacizumab treatment in combination with chemotherapy (FOLFIRI, FOLFOX or CAPIRI). Data analysis was performed in two distinct steps using nonlinear mixed-effects modeling.
In the first step of data analysis, a PK model for bevacizumab was developed and the model-predicted profiles were compared to the profiles reported by Lu et al. [26] to evaluate whether the data collected in the current study were in accordance with previous observations. Although the two-compartment model did not provide a pronounced improvement of the model fit over a one-compartment model, it was found to support the previously reported model structure and parameters and was therefore selected.
In the second step of data analysis, a TMDD model using the QSS approximation [22][23][24] was developed to characterize the in vivo bevacizumab-VEGF 165 interaction based on a simultaneous analysis of total bevacizumab and free VEGF 165 serum concentrations. The simultaneous data fitting allowed for the evaluation of the effect of bevacizumab on the reduction in free VEGF 165 levels and shed light on the in vivo binding affinity of bevacizumab to VEGF 165 as most of the current data come from in vitro biological studies [43,44]. Furthermore, it provided some valuable information on the PK properties of both bevacizumab and VEGF 165 in patients with cancer.
The TMDD model could adequately describe total bevacizumab and free VEGF 165 profiles observed during bevacizumab treatment. The simultaneous data analysis allowed for a higher precision in the estimation of PK parameters compared to the analysis of bevacizumab data alone, indicating that VEGF 165 concentrations may carry information on bevacizumab disposition. In this model, the typical volume of the central and peripheral compartment was similar to the values previously reported by Lu et al. [26] and Gaudreault et al. [45]. Bevacizumab clearance was found to be consistent with that observed by Lu et al.  [46]. The estimated free VEGF 165 levels at baseline as well as the elimination rate constant, k out , of free VEGF 165 were also in line with published values in patients with advanced cancer [47,48]. However, a much larger K ss (267 nM) than the in vitro  [1,44] was observed. This finding is in agreement with the QSS approximation that often predicts a greater K ss value in vivo than the in vitro KD [21,23].
Model-based simulations (Fig. 3) revealed a significant drop (73-80 %) in the free serum VEGF 165 concentrations upon administration of the first dose of bevacizumab, which reaches a pseudo-steady-state after multiple doses. A similar behavior was suggested by Stefanini et al. [49] in case bevacizumab is confined to the blood compartment. They mentioned a decline of 87 % in free serum VEGF 165 levels after the first dose of bevacizumab. Model predictions regarding the effect of bevacizumab on the reduction in free VEGF 165 levels in serum are of particular importance as they could indicate when dose adjustments are needed to achieve sufficient VEGF-A blockade and thus optimal anti-angiogenic activity of bevacizumab. Increases in total serum VEGF 165 concentrations were also noticed, as in previous studies by Gordon et al. [46] and Stefanini et al. [49]. This increase in total VEGF 165 could be a result of bevacizumab-VEGF 165 complex dissociation, a decrease in VEGF 165 clearance caused by the complexation process [50] or a constant production rate of VEGF 165 . Fig. 3 Model-predicted concentrations of total bevacizumab, total and free VEGF 165 for a typical patient of 70 kg. Panels a, b show the total bevacizumab and free VEGF 165 concentration profiles at doses of 5 and 7.5 mg/kg, respectively. Panel c depicts the total VEGF 165 profiles over time for the two dosing regimens The effect of demographic characteristics and genetics on the PK and VEGF-related parameters in the target population was explored. Among all covariates tested, none was found to be significant according to the predefined statistical criteria. The difficulty in identifying any significant covariate effects could be attributed to the sparse sampling schedule and the small size of the study population. Thus, the developed TMDD (binding) model did not include any covariates except for body weight in all clearance and volume parameters to ensure model stability, considering the strong biological prior [38,42] and the outcome of previous analyses [25,26]. Despite the reported gender differences in bevacizumab clearance and volume of distribution [26,45], the gender effect was not retained in the final model because gender and patient body weight are often strongly correlated, and thus, the inclusion of body weight is usually sufficient to completely describe the gender influence. The effect of serum albumin and alkaline phospatase on clearance could not be assessed in the current study since data were not available for all patients. The influence of co-administered chemotherapy (FOLFIRI, FOLFOX, CAPIRI) was not tested on bevacizumab PK, as more information would be needed to unravel potential drug-drug interactions. Nevertheless, published data indicate no significant PK interactions between bevacizumab and other antineoplastic agents [1,26].
Interestingly, a larger K ss value was observed for patients with VEGF-2578AA, VEGF-634CC and VEGF-634GC genotypes. It is noteworthy that when the covariate effect of VEGF-2578C/A or VEGF-634G/C polymorphism was added in K ss , none of the model-predicted parameters was substantially affected except for the value of K ss , which was reduced from 267 to 221 or 176 nM, respectively. This finding indicates that VEGF-2578C/A and VEGF-634G/C polymorphisms might be predictive for the binding affinity of bevacizumab to VEGF 165 , as K ss and affinity of the drug to its molecular target are correlated (Eq. 10). Larger cohort studies in bevacizumab-treated patients are needed though, to confirm these results and to assess the potential role of VEGF-2578C/A and VEGF-634G/C polymorphisms in guiding patient and dose selection.
Some limitations of our study should be also addressed. The elimination clearance of the bevacizumab-VEGF 165 complex could not be estimated and it was set equal to the clearance of the free bevacizumab, allowing for a simultaneous model fit of total bevacizumab and free VEGF 165 serum concentrations. Free bevacizumab (~150 kDa) is known to undergo proteolytic catabolism mediated by the neonatal Fc receptor (FcRn), which contributes to the slow elimination of the drug from the systemic circulation [12]. It is likely that the bevacizumab-VEGF 165 complex (assuming a 1:1 molecular interaction), which has a similar molecular weight with bevacizumab (~190 kDa) and is salvaged from degradation through binding of its Fc moiety to FcRn, exhibits a similar elimination rate. This hypothesis is further supported by Hsei et al. [50], who observed no statistically significant difference between the clearance of the free bevacizumab and the bevacizumab-VEGF 165 complex in rats. In the current study, the free drug and the complex were also assumed to be distributed in the same space, given that such large molecules would extravasate relatively slow. However, Hsei et al. [50] demonstrated a small, though significant change (20-25 %) in the volumes of distribution, implying that free bevacizumab and bevacizumab-VEGF 165 complex might not be distributed exactly in the same space. A future study in humans, which would provide information on the concentration-time profile of the bevacizumab-VEGF 165 complex, could therefore shed more light on the disposition properties of the complex and detect potential interspecies differences. Moreover, it was hypothesized that only one type of complex can be formed through monomeric binding (bevacizumab bound to one molecule of VEGF 165 ). The possibility of formation of other complexes through multimeric binding (e.g., bevacizumab bound to two molecules of VEGF 165 , VEGF 165 bound to two molecules of bevacizumab or other immune complexes) was ignored, as the bevacizumab bioassay could only detect the complex formed between bevacizumab and one molecule of VEGF 165 . The quantification of higher-order immune complexes is considered extremely difficult due to their rapid elimination from the systemic circulation [51]. Their assessment, though, could be proved valuable in reflecting an additional clearance pathway for bevacizumab. The effect of platelets on taking up VEGF and bevacizumab [52] was not included in the current model but it would be useful to be added when more data become available, as it could provide information on the fluctuation of serum VEGF levels during bevacizumab treatment.

Conclusion
In conclusion, the developed binding model adequately characterizes the PK of bevacizumab and the relationship between bevacizumab and VEGF 165 in adult patients with stage IV CRC receiving bevacizumab treatment in combination with chemotherapy. To the best of our knowledge, this is the first time a TMDD model was applied to characterize the in vivo interaction of bevacizumab with its soluble ligand, VEGF 165 . Although no significant effect of VEGF-A polymorphisms on the relationship between bevacizumab and VEGF 165 was identified, correlations between the binding affinity of bevacizumab to VEGF 165 and the VEGF-2578C/A and VEGF-634G/C SNPs were noticed. This model could serve as a basis for further studies to elucidate the role of VEGF-A polymorphisms and serum VEGF levels in treatment response and to support dose rationale for bevacizumab when combined with chemotherapy.