Exposure–response analysis using time-to-event data for bevacizumab biosimilar SB8 and the reference bevacizumab

Purpose: This analysis aimed to characterize the exposure–response relationship of bevacizumab in non-small-cell lung cancer (NSCLC) and evaluate the efficacy of SB8, a bevacizumab biosimilar, and Avastin^®, the reference bevacizumab sourced from the European Union (EU), based on the exposure reported in a comparative phase III efficacy and safety study (EudraCT, 2015-004026-34; NCT 02754882).

Materials and methods: The overall survival (OS) and progression-free survival (PFS) data from 224 patients with steady-state trough concentrations (C_ss,trough) were analyzed. A parametric time-to-event (TTE) model was developed using NONMEM^®, and the effects of treatments (SB8 and bevacizumab-EU) and patient demographic and clinical covariates on OS and PFS were evaluated. Simulations of median OS and PFS by bevacizumab C_ss,trough were conducted, and concentrations required to achieve 50% and 90% of the maximum median TTE were computed.

Results: A log-logistics model with C_ss,trough best described the OS and PFS data. Treatment was not a predictor of the hazard for OS or PFS. Simulations revealed steep exposure–response curves with a phase of rapid rise before saturating to a plateau. The median C_ss,trough values of SB8 and bevacizumab-EU reported from the clinical study were on the plateaus of the exposure–response curves. The concentrations required to achieve 50% and 90% of the maximum effect were 82.4 and 92.2 μg/mL, respectively, for OS and 79.7 and 89.1 μg/mL, respectively, for PFS.

Conclusion: Simulations based on the constructed TTE models for OS and PFS have well described the exposure–response relationship of bevacizumab in advanced NSCLC. The analysis demonstrated comparable efficacy between SB8 and bevacizumab-EU in terms of OS and PFS based on their exposure levels.

Introduction

SB8 is a biosimilar of the reference biologic bevacizumab (Avastin^®; Roche), a recombinant humanized immunoglobulin G1 monoclonal anti-vascular endothelial growth factor (VEGF) antibody. Bevacizumab inhibits VEGF-induced tumor angiogenesis and suppresses tumor growth by binding to and neutralizing the biological activity of the VEGF (Ferrara et al., 2004; Roche Registration GmbH, 2023). Bevacizumab is used in combination with chemotherapy to treat various solid tumor types. Its therapeutic indications vary slightly across countries. In the European Union (EU), bevacizumab is authorized for the treatment of metastatic colorectal cancer, metastatic breast cancer, non-small-cell lung cancer (NSCLC), metastatic renal cell carcinoma, epithelial ovarian cancer, fallopian tube cancer, primary peritoneal cancer, and cervical cancer (Roche Registration GmbH, 2023).

SB8 was approved in the EU in 2020/2021 (Aybintio^®/Onbevzi™) for the same cancer types as for Avastin^®, based on the demonstration of comparable structural, functional, non-clinical, clinical pharmacokinetic (PK), efficacy, and safety profiles to the reference bevacizumab. The product subsequently received approval in the Republic of Korea (Onbevzi™), in Canada (Aybintio^®) in 2021, and in Taiwan (Onbevzi™) in 2022 (Taiwan Food and Drug Administration, 2018; IPPR, 2021; Samsung Bioepis NL B.V, 2023b; Samsung Bioepis NL B.V, 2023a; Samsung Bioepis Co., Ltd., 2023). Analytical testing showed similarity between SB8 and the reference bevacizumab in terms of physicochemical properties (e.g., structure, biological activity, purity, and stability), and non-clinical testing showed functional similarity and in vivo pharmacological/toxicological activity between the two drug products. In support of establishing clinical similarity, two clinical studies were undertaken to compare the PK, efficacy, and safety of SB8 with its reference bevacizumab. In a phase I, randomized, single-dose, three-arm comparative PK study of healthy volunteers, the PK similarity of SB8 to both EU-sourced bevacizumab (bevacizumab-EU) and United States (US)-sourced bevacizumab (bevacizumab-US) was demonstrated, with the 90% confidence interval (CI) of the test-to-reference geometric mean ratio for the area under the serum drug concentration–time curves, as well as the maximum concentration, falling between 80% and 125% (Shin et al., 2020). In addition to the comparative phase I study, a phase III, randomized, double-blind, multicenter comparative study was conducted to evaluate the efficacy, safety, PK, and immunogenicity of SB8 compared with bevacizumab-EU in patients with metastatic or recurrent NSCLC. The study showed equivalence in the best overall response rate risk ratio, with comparable safety, PK, and immunogenicity profiles between SB8 and bevacizumab-EU (Reck et al., 2020). Upon the product’s approval, network meta-analysis, which combines direct and indirect evidence for the simultaneous comparison of multiple medicinal products, has been conducted, showing the similarity of SB8 to the reference product and other bevacizumab biosimilars (Xu et al., 2022).

The nature and scope of these studies performed during the development of SB8 offer a glimpse into the regulatory framework of biosimilar products, which is quite different from that associated with the development pathway for novel biologic drugs. The focus of the biosimilar development program is not to establish patient benefit per se but to demonstrate comparability at the analytical, functional, and clinical levels, as the development of biosimilars relies on existing scientific knowledge about the safety and effectiveness of the approved reference products (Weise et al., 2012). Nonetheless, conceptually, PK/PD (pharmacodynamic) modeling and simulation during biosimilar development can play a role similar to its role in the development of novel therapeutic proteins (Dodds et al., 2013; Wang et al., 2019). Pharmacometric analysis can provide insights into the exposure range that would be sensitive enough to assess clinically meaningful differences between biosimilars and reference products and can optimize study designs accordingly. However, such a similarity approach should be preceded by an understanding of the well-established exposure–response relationship of reference biologics and any variability associated with treatment.

As with most targeted protein therapeutics in oncology, patient responses may vary following treatment with bevacizumab in combination with chemotherapy. Multiple factors can contribute to this inter-patient variability, including, but not limited to, PK/biodistribution, drug interactions, tumor burden, and tumor heterogeneity (Symeonides and Gourley, 2015; Chakrabarti and Michor, 2017; Qi and Cao, 2023). The total plasma or serum VEGF concentration is used as an in vivo marker, and to date, its true association with the efficacy of bevacizumab for treating cancer in humans is uncertain, as the results of studies assessing the predictive value of serum VEGF levels for treatment outcomes have been not consistent (Presta et al., 1997; Bernaards et al., 2010; Tang et al., 2016). Aiming to find biomarkers, several studies have investigated the exposure–response relationship of bevacizumab, exploring trough concentrations associated with a threshold for treatment efficacy (Nugue et al., 2013; Caulet et al., 2016; Akbulut et al., 2018; Papachristos et al., 2020), and most of the research has focused on metastatic colorectal cancer, with little to no published data from investigations of bevacizumab’s exposure–response relationship in the treatment of advanced NSCLC.

Therefore, our modeling and simulation focused on the exposure–response modeling of bevacizumab, specifically for its use in the treatment of metastatic or recurrent NSCLC. This study aimed to 1) characterize the exposure–response relationship of SB8 and bevacizumab-EU in patients with advanced NSCLC and 2) evaluate the therapeutic efficacy of SB8 and bevacizumab-EU based on the exposure of each drug product.

Materials and methods

Study population

The clinical data used in this analysis were collected from 224 patients in a phase III comparative efficacy and safety study (SB8-G31-NSCLC) comparing SB8 (bevacizumab biosimilar) and the reference bevacizumab in the treatment of metastatic or recurrent non-squamous NSCLC (Reck et al., 2020). The patients in this study were randomized to receive SB8 or the reference bevacizumab (Avastin^®, bevacizumab-EU), each at 15 mg/kg intravenously (IV), with paclitaxel 200 mg/m² and carboplatin area under the curve 6, every 3 weeks for 4–6 cycles of the induction period. The tumor size was assessed by radiographic imaging, and patients with complete response, partial response, or stable disease as per the Response Evaluation Criteria in Solid Tumors (RECIST) after the induction period were given maintenance therapy with SB8 or bevacizumab-EU, each at 15 mg/kg IV every 3 weeks until progressive disease, unacceptable toxicity, death, or the end of the study, whichever occurred first. Trough serum concentrations of SB8 and bevacizumab-EU were assessed by pre-dose sampling before IV infusion at cycles 1, 3, 5, and 7. The exposure–response analysis included the overall survival (OS) and progression-free survival (PFS) assessed by blinded independent central review. Data from patients with available steady-state trough concentrations for SB8 (n = 100) and bevacizumab-EU (n = 124) measured just before cycle 5 or cycle 7 were included in the analysis.

Study SB8-G31-NSCLC was registered at EudraCT (2015-004026-34) and ClinicalTrials.gov (NCT 02754882). The study was in compliance with the Declaration of Helsinki, International Conference on Harmonization (ICH), and Good Clinical Practice (GCP) guidelines, as well as applicable local regulations. The protocol was approved by the local ethics committees at each study center, and informed consent was obtained from all patients before enrolling in the trial.

Parametric TTE model development

Modeling efforts were streamlined using graphical exploratory analyses of OS and PFS with the Kaplan–Meier method before building a base model. A non-parametric log-rank test was performed for the unadjusted OS and PFS analyses according to the treatment administered. Cox proportional hazards regression modeling, with treatment (SB8 and bevacizumab-EU) as an explanatory variable, was used to estimate hazard ratios and associated 95% CIs.

OS and PFS were described using a TTE modeling approach (Lim, 2021). Base models for OS and PFS without any covariates were created and selected first, and potential covariates were assessed thereafter. Various TTE models were tested, including exponential (Eq. 1), Weibull (Eq. 2), and log-logistic (Eq. 3) distributions of event times, and a proportional hazard model was used for the covariate effect. The final base models for OS and PFS were selected based on the lowest objective function value (OFV), graphical and statistical evaluations, and scientific plausibility.

${S\left(t\right)=e}^{-\lambda t},(1)$

$S\left(t\right)=e^{-{\left(\lambda t\right)}^{\gamma }},(2)$

$S\left(t\right){=\left\{1+{\left(\frac{t}{\lambda }\right)}^{\gamma }\right\}}^{-1},(3)$

where S(t) is a survival function, $\lambda$ is the degree of decrease in S(t) over time, and $\gamma$ is an additional parameter providing more flexibility in S(t) over time.

The goodness of fit of the models to the observed data was evaluated using statistical methods and the visual predictive check (VPC). The Kaplan–Meier plots of the observed data were summarized as median estimates and 95% CIs, and then, they were compared with the simulated medians and 95% prediction intervals (PIs) (2.5th–97.5th percentiles) determined using the respective models for OS and PFS. The VPC values were further stratified by quartiles of the observed bevacizumab C_ss,trough. The Wald test was used for model evaluations, and the likelihood ratio test was used for comparing hierarchical models. A decrease of at least 3.84, corresponding to a p-value of 0.05, was considered a statistically significant improvement of fit during base model construction.

Covariate modeling was performed by the forward stepwise addition of covariates to the model, followed by a backward elimination procedure from the full model. In the univariate analysis, the statistical significance of each covariate was tested individually with a criterion for changes in an OFV of 6.63 (α = 0.01; df = 1). Covariates with the most significant decreases in the OFV were collectively set as a new base covariate model, and the process was repeated until there were no further covariates indicating a significant decrease in the OFV. If two or more covariates were included upon attainment of the full model, a more stringent OFV criterion of 7.88 (α = 0.005, df = 1) was applied to allow covariates to be maintained in the model.

The following covariates describing patient characteristics were selected for screening based on physiological plausibility and clinical relevance.

Age (<65 years vs. ≥65 years); sex (male vs. female); race (white vs. non-white); body weight (continuous); body mass index (BMI) (continuous); subtype of lung cancer (adenocarcinoma, large-cell neuroendocrine carcinoma, large-cell carcinoma, adenosquamous carcinoma, or other); Eastern Cooperative Oncology Group (ECOG) performance status at enrollment (0 vs. 1); smoking history (never smoker vs. former smoker vs. current smoker); treatment (SB8 vs. bevacizumab-EU); overall ADA (anti-drug antibody) results up to cycle 7 and overall ADA results by the end of treatment (positive vs. negative vs. inconclusive); pre-dose steady-state (continuous) concentrations of bevacizumab (SB8 and bevacizumab-EU); geographical region of residence (EU vs. non-EU); and country (Russia vs. Ukraine vs. Georgia; Hungary; Germany; Spain vs. Belarus; Romania; Serbia; Poland vs. Korea; Thailand; and Taiwan).

The demographic and clinical characteristics of the patients are summarized in Table 1.

TABLE 1

TABLE 1. Demographic and clinical characteristics of patients included in the time-to-event analysis.

A proportional hazard model was implemented to test each of the covariates on the baseline hazard. C_ss,trough was taken from the pre-dose serum concentrations of SB8 and bevacizumab-EU measured mostly before cycle 7 (221 samples) or else cycle 5 (3 samples) at ≥100 days after the first dosing in each patient (Roche Registration GmbH, 2005). To describe the drug’s concentration effect on survival or disease progression in terms of fractional decrease of the hazard, an inhibitory sigmoidal maximum effect (E_max) model (Eq. 4) was implemented as follows (Chigusta et al., 2017):

${H\left(t\right)=H}_0\left(t\right)\left[1-\left(\frac{E_{\mathit{max}}*{C_{ss,trough}}^{Hill}}{{{EC}_{50}}^{Hill}+{C_{ss,trough}}^{Hill}}\right)\right],(4)$

where H(t) and ${\mathrm{H}}_0$ (t) are the hazard and baseline hazard, respectively; $C_{ss,trough}$ is the steady-state trough concentration of treatment; $E_{\mathit{max}}$ is the maximum inhibitory effect of treatment on the hazard; ${EC}_{50}$ is $C_{ss,trough}$ at half of $E_{\mathit{max}}$; and Hill is the hill coefficient.

For categorical covariates with missing data, a mixture model was applied in NONMEM^®, version 7.4 (ICON Development Solutions, Ellicott City, MD, United States), to impute missing data based on weighted maximum likelihood (Frame et al., 2003). In the mixture model, the missing proportion was fixed at the respective proportion of the observed data under the assumption that the data were missing completely at random.

Simulation study for the exposure–response relationship

The exposure–response curve for OS and PFS by the range of C_ss,trough was simulated by implementing inverse functions for the TTE functions in NONMEM codes for the respective final models.

To further determine concentrations corresponding to 50% of the simulated maximum median time to event, where the event is death for OS and disease progression or death for PFS (hereafter abbreviated as TTE50CONC), 50% of the maximum median TTE was first calculated from the half of the difference between the simulated baseline median TTE and the maximum median TTE (plus the baseline median TTE value) estimated using the final models for OS and PFS. The median TTE associated with concentrations ranging from 0 to 150 μg/mL was tested at increments of 0.1 μg/mL during simulation. The conditional statement in the $ERROR block in the NM-TRAN control records allowed the tested concentration to be capped if the response associated with that concentration reached the half-maximal response in the exposure–response curve (Supplementary Table S1). This way, the largest concentration approaching the half-maximal response was determined to be TTE50CONC. Similarly, concentrations corresponding to 90% of the simulated maximum median TTE for OS and PFS (hereafter abbreviated as TTE90CONC) were derived using the same approach to assess 90% of the maximal response.

An overview of the workflow for this modeling and simulation is provided in Figure 1.

FIGURE 1

FIGURE 1. Workflow for time-to-event simulation modeling and the analysis of overall survival and progression-free survival. OS, overall survival; PFS, progression-free survival.

All model development and parameter estimations throughout this study were performed by first-order conditional estimation with the interaction method using NONMEM. Data handling, graphical model diagnosis, and statistical analyses were performed using R, version 4.2.3 (R Foundation for Statistical Computing, Vienna, Austria).

Results

The exploratory analyses for OS and PFS using the combined data (SB8 + bevacizumab-EU) and separate data by the treatment group show comparable survival curves between SB8 and bevacizumab-EU (Figure 2). The log-rank test revealed no treatment differences in terms of both OS (p = 0.6) and PFS (p = 0.6), and the estimated hazard ratios were 1.13 (95% CI: 0.74–1.74; p = 0.6) for OS and 1.07 (95% CI: 0.78–1.48; p = 0.7) for PFS.

FIGURE 2

FIGURE 2. Kaplan–Meier plot for (A) overall survival for SB8 + BEV-EU; (B) overall survival by the treatment group; (C) progression-free survival for SB8 + BEV-EU; and (D) progression-free survival by Cox proportional hazard model results for OS and PFS analyzed by SB8 and bevacizumab-EU treatment groups. BEV, reference EU-sourced bevacizumab; HR, hazard ratio; OS, overall survival; and PFS, progression-free survival.

The log-logistic distribution model best described the time to OS and PFS. Based on the log-likelihood ratio test, the Weibull distribution model was superior to the exponential distribution model in all cases. The log-logistic distribution model had the lowest OFV for both OS and PFS. Although testing for statistical significance between Weibull and log-logistic was not applicable since these models do not have a full-nested relationship, the predictive performance for OS and PFS was similar between the log-logistic and Weibull models, as evidenced by VPC plots. The VPC results show the 95% PI of the simulated data overlaid on the Kaplan–Meier curve of the observed data, indicating that the model predictions aligned well with the observed data (Figure 3). The median, 95% PI, and 95% CI largely overlapped each other on VPC results for OS and PFS stratified by quartiles of bevacizumab C_ss,trough (Figure 4, Figure 5).

FIGURE 3

FIGURE 3. Visual predictive check plots for final models of overall survival (top) and progression-free survival (bottom). CI, confidence interval; KM, Kaplan–Meier; and PI, prediction interval.

FIGURE 4

FIGURE 4. Visual predictive check plots for overall survival stratified by exposure quartiles: (A) C_ss,trough quartile 1; (B) C_{ss trough} quartile 2; (C) C_ss,trough quartile 3; and (D) C_ss,trough quartile 4. CI, confidence interval; KM, Kaplan–Meier; and PI, prediction interval.

FIGURE 5

FIGURE 5. Visual predictive check plots for progression-free survival stratified by exposure quartiles: (A) C_ss,trough quartile 1; (B) C_ss,trough quartile 2; (C) C_ss,trough quartile 3; and (D) C_ss,trough quartile 4. CI, confidence interval; KM, Kaplan–Meier; and PI, prediction interval.

C_ss,trough was a significant predictor of the hazard for OS and PFS. Covariate analysis did not indicate any other patient demographic or clinical characteristics to be significant enough to be included in the final model (Table 2). In the univariate analysis for the OS model, BMI, weight, and C_ss,_trough were significant based on the criterion of p = 0.01 (ΔMOFV = −6.63), but BMI and weight did not meet the criteria during the forward selection procedure and, thus, were not included in the full model. In the univariate analysis for the PFS model, only C_ss,trough was significant based on the criterion of p = 0.01 (ΔMOFV = −6.63) and was, thus, included in the final model. Treatment (SB8 vs. bevacizumab-EU) did not affect the models for OS and PFS. C_ss,trough was not a significant covariate for any model when a proportional hazard model was used, but when a sigmoid maximum effect model was implemented, C_ss,trough significantly improved the model. The Hill coefficient in the sigmoid maximum effect model was tested with a range of 1–100, and a value of 20 was fixed in the model. The selected value of the Hill coefficient suggests a nearly all-or-none effect of SB8 and bevacizumab-EU on OS and PFS. The C_ss,trough values at half of E_max (EC₅₀) for the OS and PFS models were estimated as 77.5 and 77.7 μg/mL, respectively. The final TTE models for both OS and PFS are described in Table 3.

TABLE 2

TABLE 2. Model selection for the time-to-event analysis.

TABLE 3

TABLE 3. Parameter estimates in the final models for overall survival and progression-free survival.

Simulations of median OS and PFS by ranges of C_ss,trough were carried out using the constructed final TTE models. A steep exposure–response relationship was shown, wherein the predicted median OS and PFS rise rapidly with increasing C_ss,trough, plateau after a certain value (i.e., the inflection point), and then, remain unchanged with further increases in C_ss,trough (Figure 6). When comparing the interquartile range of C_ss,trough to the exposure–response curve, the lower end of the range for both treatments was on the rising curve, while the upper range of concentrations was at a plateau. The median C_ss,trough values of SB8 (115.114 μg/mL) and bevacizumab-EU (125.453 μg/mL) were observed from the clinical study of the simulations; both treatments were well above the steep part of the exposure–response curves, and the 95% bootstrap median CIs of C_ss,trough for SB8 and bevacizumab-EU overlapped one another on the plateaus. The 50th and 90th percentiles of the predicted median OS were 17.1 months and 21.2 months, respectively. The 50th and 90th percentiles of the predicted median PFS were 9.38 months and 10.3 months, respectively. Based on these values, the TTE50CONC and TTE90CONC values derived from simulation were 82.4 and 92.2 μg/mL, respectively, for OS and 79.7 and 89.1 μg/mL, respectively, for PFS.

FIGURE 6

FIGURE 6. Exposure–response curves for overall survival and progression-free survival based on the final models. The vertical solid lines indicate observed median steady-state trough concentrations (C_ss,trough) for SB8 (blue) and bevacizumab-EU (red) at 15 mg/kg intravenous infusion every 3 weeks. The shaded areas indicate the 95% median CI of steady-state trough concentrations for SB8 (light blue) and bevacizumab-EU (pink).

Discussion

This simulation study provided a unique modeling and simulation framework for characterizing the relationship between bevacizumab C_ss,trough and its efficacy in terms of OS and PFS in advanced NSCLC. To our knowledge, this is the first study that proceeded to publication after the conduct of a TTE model-based simulation approach and evaluation of the efficacy of a bevacizumab biosimilar along with a reference biologic based on the concentration–response (time to death and time to disease progression) relationship.

TTE modeling was performed using OS and PFS data for patients with metastatic or recurrent non-squamous NSCLC who had available pre-dose steady-state concentration data from a comparative clinical efficacy and safety study. Covariate analysis evaluated the impact of key patient-level factors on survival or disease progression. The final models for both OS and PFS were selected, with C_ss,trough as a covariate showing sigmoidal maximum effects. Simulation to assess exposure in relation to median OS and PFS was carried out by implementing the inverse functions for median TTE functions in the respective final models. The bevacizumab exposure–response curve for OS and PFS based on the corresponding final model reflects a steep relationship where the median TTE (OS or PFS) corresponds with a phase of sharp increase, resulting in a wide range of bevacizumab concentrations showing maximum efficacy.

Clinical factors affecting the survival of patients with advanced NSCLC are worth mentioning. Through the univariate analysis of covariates, we noted that the hazard for OS increased with decreasing BMI and weight; however, this increase was not statistically significant. This observation was consistent with previous findings of low BMI and significant weight loss (≥5%) as related factors associated with poor survival among patients with NSCLC (Hoang et al., 2012; Shepshelovich et al., 2019). Weight-based bevacizumab dosing for the treatment of solid cancer types, accounting for the body weight’s contribution to clearance and volume of distribution, may explain these results (Lu et al., 2008; Roche Registration GmbH, 2023). Nevertheless, the association between BMI and survival outcomes with bevacizumab treatment (along with other patient factors interacting with BMI affecting poor prognosis) remains to be fully elucidated. Furthermore, the impact of overall positive ADA results on OS and PFS showed no significance in our model, and in the comparative phase III study, ADA formation at each time point of the cycle was comparable between SB8 and bevacizumab-EU (Samsung Bioepis NL B.V, 2023b). It is to be noted that the clinical significance of anti-bevacizumab is not known (FDA, 2022). In neither of the covariate analyses for OS and PFS did the treatment (SB8 or bevacizumab-EU) affect hazards differently, suggesting that the treatment effects on OS and PFS were comparable between SB8 and bevacizumab-EU.

In the final model, an E_max model was used to empirically describe hazard reduction according to a range of concentrations (C_ss,trough). While this model may be limited in its physiological basis, the relationship between the concentration of a drug and its pharmacological effects, whether the variable is binary or continuous, can be conveniently described by a sigmoid E_max model, especially in cases where the relationship is close to that of a cumulative log-normal distribution (Bowling et al., 2009; Proost et al., 2020). In finding the best-fit model, a wide range of Hill coefficients was tested, from a simple E_max model (γ = 1) to a steep model (γ = 100), to describe the drug’s pharmacological effect. This way, a proper estimation of parameters concerning pharmacological efficacy (E_max) and potency (EC₅₀) was obtained to determine the hazard function according to drug concentrations. When a hazard was determined based on a given concentration value in the final TTE model, the median OS and PFS were subsequently determined from the corresponding hazard, and TTE50CONC and TTE90CONC were further determined from simulation to provide clinical context to the hazard.

As the last step in model development, the TTE VPC was used to evaluate model performance, showing that the final models for OS and PFS fit well with the observed survival data in the Kaplan–Meier curves. However, around year 2, the model predicted a slightly lower median survival value than the observed median values. Nonetheless, the 95% CI for the observed values overlapped with the 95% PI throughout the observation period, indicating no statistical significance. Also, the low number of patients in the data at this time point influenced diagnostic accuracy and precision.

The simulated exposure–response curve aligns well with the response associated with steady-state concentrations measured following multiple doses of bevacizumab at 7.5 and 15 mg/kg, as reported in historical investigations of the reference product, Avastin^® (FDA, 2004; Johnson et al., 2004; FDA, 2006). With C_ss,trough between 73 ± 43 μg/mL and 135 ± 48 μg/mL reported for the respective doses of 7.5 and 15 mg/kg of Avastin^® in a phase II NSCLC study (AVF0757g), an optimal therapeutic response may be seen in many patients with a dose of 7.5 mg/kg; concentrations within this range are right at or above those located at the rapid increase of the median TTE curve. Moreover, when concentrations are increased with a higher dose at 15 mg/kg, the response associated with the upper end of the reported concentration range (i.e., 135 μg/mL) would show the maximum efficacy of bevacizumab based on the exposure–response curves. Similarly, the interquartile range of C_ss,trough (92.0–150.5 μg/mL) for SB8 and bevacizumab-EU was associated with the rising and plateau of the exposure–response curve, indicating therapeutic response. These comparisons show that the modeling and simulation accurately and predictably characterized the concentration-dependent response of bevacizumab within the therapeutic range.

In addition to visual representations of the exposure–response relationship, the shape of the curve was further explored to understand the range of relevant values by computing the concentrations corresponding to 50% and 90% of the maximum median TTE (TTE50CONC, TTE90CONC) for OS and PFS, respectively. The median and interquartile range of median C_ss,trough of both SB8 and bevacizumab-EU, along with the previously reported range for C_ss,trough of Avastin^®, were all well above TTE50CONC. While other methods, such as the use of inverse functions of the simulation codes in the final model, may be feasible for estimation, this requires the analytical calculation of each relevant parameter in the TTE function and Hill equations to derive the final value. Thus, in our analyses, conditional statements were used to find concentrations associated with the response of interest and then to cap the values once they were determined. For this purpose, the simulation dataset was prepared with a narrow concentration interval ranging from zero to the maximum bevacizumab C_ss,trough reported in the clinical study, in increments of 0.1 μg/mL, enabling NONMEM to compute TTE50CONC and TTE90CONC with reasonable accuracy.

Our modeling analyses based on the principle of TTE have utility in evaluating the comparative efficacy of two treatments. The analyses remain empirical and exploratory for the purpose of predicting the effect of dose changes and dose optimization, as the model is based on clinical data for monoclonal antibodies tested at a single-dose level with available C_ss,trough. Further improvement can be made by applying a mixed-effects analytical approach. First, the effects of inter-individual variabilities of bevacizumab PK were not incorporated, and thus, typical prediction of TTE outcomes was made from the current final model. Estimating and incorporating the inter-individual variability of bevacizumab PK into TTE models help elicit valuable insights into a predictive TTE response interval in relation to drug exposure, which in turn provides more biologically plausible predictions for OS and PFS based on individual plasma concentration levels. Also, studies using PK/PD modeling may discern additional concentration-dependent covariates and time-varying PK parameters (e.g., time-varying clearance with tumor dynamics), potentially affecting the final TTE models for OS and PFS (Liu et al., 2017). Thus, given the sensitivity of tumor size to monoclonal antibodies used in oncology and the association between tumor growth rate and survival in the context of advanced NSCLC (Claret et al., 2009; Wang et al., 2009; Boyle, 2010; Gong et al., 2020), it would be worth assessing the impact of changes in tumor size after bevacizumab treatment. In particular, a separate model for predicting the tumor size can be developed to integrate model-predicted changes in the tumor size into the survival model, and the current model may be updated based on early readouts related to tumor size. Ultimately, such analysis based on tumor size-integrated TTE modeling, along with modeling and simulation results such as those generated by our study, could guide in predicting outcomes associated with the same classes of anti-VEGF drugs using early readouts on tumor size.

From the perspective of biosimilarity evaluation, our study extended the application of PK/PD modeling and simulation to evaluate the comparative efficacy of a biosimilar and its reference biologic based on their exposures (Wang et al., 2019). In a biosimilar development program, modeling and simulation are typically applied to optimize study design, ensuring adequate sensitivity in the selected study population, appropriate dosing, and sufficient sample size based on an understanding of inter-individual variability in PK/PD (Zhu et al., 2018; Li et al., 2020). In addition, a quantitative relationship between a PD biomarker relevant to a drug’s mechanism of action, as well as clinical endpoints, can be assessed through the PK/PD model, providing the information necessary to assess the model’s utility as a surrogate primary endpoint in clinical research (Lee et al., 2011; Li L et al., 2018; Jang et al., 2022). In cases where limited comparative data are available for biosimilars, network meta-analysis can be performed to assess the similarity of the products (Chiumente et al., 2017), and such analysis may be conducted through modeling as well. In addition to these well-known, conventional applications, our exposure–response analysis using a TTE model provides a novel “science-applied tool” for assessing the potential impact of exposure differences between a biosimilar and the reference product on clinical efficacy.

In conclusion, our modeling and simulation analysis elucidated the exposure–response relationship of bevacizumab in the treatment of patients with advanced NSCLC, providing fundamental information about the drug. This relationship further substantiates the similarity in treatment effects between SB8 and EU-bevacizumab.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

SP: writing–original draft, writing–review and editing, formal analysis, methodology, project administration, and validation. JJ: writing–review and editing. SJ: writing–review and editing. H-SL: writing–review and editing, conceptualization, formal analysis, methodology, and validation.

Funding

The authors declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

The authors thank Samsung Bioepis Co., Ltd., for providing the clinical data used for this research.

Conflict of interest

Authors JJ and SJ were employed by Samsung Bioepis Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1288308/full#supplementary-material

References

Akbulut, H., Ocal, M., Sonugur, F. G., Abgarmi, S. A., Babahan, C., Akay, B., et al. (2018). The trough levels of bevacizumab significantly affect the outcome of the treatment in patients with metastatic colorectal cancer: a Turkish Oncology Group study. Clin. Oncol. 36 (15), e15553. doi:10.1200/JCO/2018.36.15_suppl.e15553

CrossRef Full Text | Google Scholar

Bernaards, C., Hegde, P., Chen, D., Holmgren, E., Zheng, M., Koeppen, A. M., et al. (2010). Circulating vascular endothelial growth factor (VEGF) as a biomarker for bevacizumab-based therapy in metastatic colorectal, non-small cell lung, and renal cell cancers: analysis of phase III studies. Clin. Oncol. 28 (15), 10519. doi:10.1200/jco.2010.28.15_suppl.10519

CrossRef Full Text | Google Scholar

Bowling, S. R., Khasawneh, M. T., Sittichai, K., and Cho, B. R. (2009). A logistic approximation to the cumulative normal distribution. J. Ind. Eng. Manag. 2 (1), 114–127. doi:10.3926/jiem.2009.v2n1.p114-127

CrossRef Full Text | Google Scholar

Boyle, R. M. (2010). The use of biologics in cancer therapy. U. S. Pharm. 35 (3), 4–7. Avaliable at: https://uspharmacist.com/article/the-use-of-biologics-in-cancer-therapy (Accessed August 12, 2023).

Google Scholar

Caulet, M., Lecomte, T., Bouché, O., Rollin, J., Gouilleux-Gruart, V., Azzopardi, N., et al. (2016). Bevacizumab pharmacokinetics influence overall and progression-free survival in metastatic colorectal cancer patients. Clin. Pharmacokinet. 55 (11), 1381–1394. doi:10.1007/s40262-016-0406-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakrabarti, S., and Michor, F. (2017). Pharmacokinetics and drug interactions determine optimum combination strategies in computational models of cancer evolution. Cancer Res. 77 (14), 3908–3921. doi:10.1158/0008-5472.CAN-16-2871

PubMed Abstract | CrossRef Full Text | Google Scholar

Chigutsa, E., Long, A. J., and Wallin, J. E. (2017). Exposure-response analysis of necitumumab efficacy in squamous non-small cell lung cancer patients. CPT Pharmacometrics Syst. Pharmacol. 6 (8), 560–568. doi:10.1002/psp4.12209

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiumente, M., Mengato, D., and Messori, A. (2017). Rituximab biosimilar evaluated by network meta-analysis. Haematologica 102 (12), e497–e498. doi:10.3324/haematol.2017.175042

PubMed Abstract | CrossRef Full Text | Google Scholar

Claret, L., Girard, P., Hoff, P. M., Van Cutsem, E., Zuideveld, K. P., Jorga, K., et al. (2009). Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics. J. Clin. Oncol. 27 (25), 4103–4108. doi:10.1200/JCO.2008.21.0807

PubMed Abstract | CrossRef Full Text | Google Scholar

Dodds, M., Chow, V., Markus, R., Perez-Ruixo, J. J., Shen, D., and Gibbs, M. (2013). The use of pharmacometrics to optimize biosimilar development. J. Pharm. Sci. 102 (11), 3908–3914. doi:10.1002/jps.23697

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrara, N., Hillan, K. J., Gerber, H. P., and Novotny, W. (2004). Discovery and development of bevacizumab, an anti-VEGF antibody for treating cancer. Nat. Rev. Drug Discov. 3, 391–400. doi:10.1038/nrd1381

PubMed Abstract | CrossRef Full Text | Google Scholar

Food and Drug Administration (2004). Bevacizumab clinical Pharmacology and biopharmaceutical review. Available at: https://www.accessdata.fda.gov/drugsatfda_docs/nda/2004/STN-125085_Avastin_BioPharmr.pdf (Accessed August 13, 2023).

Google Scholar

Food and Drug Administration (2006). Bevacizumab approval package for application number: BLA 125085/85. Available at: https://www.accessdata.fda.gov/drugsatfda_docs/bla/2006/125085_0085_avastin.pdf (Accessed December 15, 2023)

Google Scholar

Food and Drug Administration (2022). Bevacizumab prescribing information. Available at: https://www.accessdata.fda.gov/drugsatfda_docs/label/2022/125085s340lbl.pdf (Accessed December 15, 2023)

Google Scholar

Frame, B., Miller, R., and Lalonde, R. L. (2003). Evaluation of mixture modeling with count data using NONMEM. J. Pharmacokinet. Pharmacodyn. 30 (3), 167–183. doi:10.1023/a:1025564409649

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, Y., Mason, J., Shen, Y. L., Chang, E., Kazandijian, D., Blumenthal, G. M., et al. (2020). An FDA analysis of the association of tumor growth rate and overall and progression-free survival in metastatic non-small cell lung cancer (NSCLC) patients. J. Clin. Oncol. 38 (15), 9541. doi:10.1200/JCO/2020.38.15_suppl.9541

CrossRef Full Text | Google Scholar

Hoang, T., Dahlberg, S. E., Sandler, A. B., Brahmer, J. R., Schiller, J. H., and Johnson, D. H. (2012). Prognostic models to predict survival in non-small-cell lung cancer patients treated with first-line paclitaxel and carboplatin with or without bevacizumab. J. Thorac. Onol 7 (9), 1361–1368. doi:10.1097/JTO.0b013e318260e106

PubMed Abstract | CrossRef Full Text | Google Scholar

International Pharmaceutical Regulators Programme (IPPR) (2021). Onbevzi public assessment summary information for biosimilar IPRP biosimilars working group. Available at: https://www.mfds.go.kr/com/file/down.do?dnCd=eng&fileNm=PASIB_Onbevzi(MFDS).pdf&filePath=contents/ (Accessed August 12, 2023).

Google Scholar

Jang, J. H., Gomez, R. D., Bumbea, H., Nogaieva, L., Wong, L. L. L., Lim, S. M., et al. (2022). A phase III, randomised, double-blind, multi-national clinical trial comparing SB12 (proposed eculizumab biosimilar) and reference eculizumab in patients with paroxysmal nocturnal haemoglobinuria. EJHaem 4 (1), 26–36. doi:10.1002/jha2.632

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D. H., Fehrenbacher, L., Novotny, W. F., Herbst, R. S., Nemunaitis, J. J., Jablons, D. M., et al. (2004). Randomized phase II trial comparing bevacizumab plus carboplatin and paclitaxel with carboplatin and paclitaxel alone in previously untreated locally advanced or metastatic non-small-cell lung cancer. J. Clin. Oncol. 22 (11), 2184–2191. doi:10.1200/JCO.2004.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. W., Jang, J. H., Kim, J. S., Yoon, S. S., Lee, J. H., Kim, Y. K., et al. (2011). Uncontrolled complement activation and the resulting chronic hemolysis as measured by LDH serum level at diagnosis as predictor of thrombotic complications and mortality in a large cohort of patients with paroxysmal nocturnal hemoglobinuria (PNH). Blood 118 (21), 3166. doi:10.1182/blood.V118.21.3166.3166

CrossRef Full Text | Google Scholar

Li, C. S. W., Sweeney, K., and Cronenberger, C. (2020). Population pharmacokinetic modeling of PF-06439535 (a bevacizumab biosimilar) and reference bevacizumab (Avastin^®) in patients with advanced non-squamous non-small cell lung cancer. Cancer Chemother. Pharmacol. 85 (3), 487–499. doi:10.1007/s00280-019-03946-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Ma, L., Schrieber, S. J., Rahman, N. A., Deisseroth, A., Farrell, A. T., et al. (2018). Quantitative relationship between AUEC of absolute neutrophil count and duration of severe neutropenia for G-CSF in breast cancer patients. Clin. Pharmacol. Ther. 104 (4), 742–748. doi:10.1002/cpt.991

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, H. S. (2021). Brief introduction to parametric time to event model. Transl. Clin. Pharmacol. 29 (1), 1–5. doi:10.12793/tcp.2021.29.e7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Yu, J., Li, H., Liu, J., Xu, Y., Song, P., et al. (2017). Association of time-varying clearance of nivolumab with disease dynamics and its implications on exposure response analysis. Clin. Pharmacol. Ther. 101 (5), 657–666. doi:10.1002/cpt.656

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, J. F., Bruno, R., Eppler, S., Novotny, W., Lum, B., and Gaudreault, J. (2008). Clinical pharmacokinetics of bevacizumab in patients with solid tumors. Cancer Chemother. Pharmacol. 62 (5), 779–786. doi:10.1007/s00280-007-0664-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Nugue, G., Bidart, M., Arlotto, M., Mousseau, M., Berger, F., and Pelletier, L. (2013). Monitoring monoclonal antibody delivery in oncology: the example of bevacizumab. PloS one 8 (8), e72021. doi:10.1371/journal.pone.0072021

PubMed Abstract | CrossRef Full Text | Google Scholar

Papachristos, A., Kemos, P., Kalofonos, H., and Sivolapenko, G. (2020). Correlation between bevacizumab exposure and survival in patients with metastatic colorectal cancer. Oncologist 25 (10), 853–858. doi:10.1634/theoncologist.2019-0835

PubMed Abstract | CrossRef Full Text | Google Scholar

Presta, L. G., Chen, H., O’Connor, S. J., Chisholm, V., Meng, Y. G., Krummen, L., et al. (1997). Humanization of an anti-vascular endothelial growth factor monoclonal antibody for the therapy of solid tumors and other disorders. Cancer Res. 57 (20), 4593–4599.

PubMed Abstract | Google Scholar

Proost, J. H., Eleveld, D. J., and Struys, MMRF (2020). Population pharmacodynamic modeling using the sigmoid Emax model: influence of inter-individual variability on the steepness of the concentration-effect relationship. A simulation study. AAPS J. 23 (1), 10. doi:10.1208/s12248-020-00549-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, T., and Cao, Y. (2023). Dissecting sources of variability in patient response to targeted therapy: anti-HER2 therapies as a case study. Eur. J. Pharm. Sci. 186, 106467. doi:10.1016/j.ejps.2023.106467

PubMed Abstract | CrossRef Full Text | Google Scholar

Reck, M., Luft, A., Bondarenko, I., Shevnia, S., Trukhin, D., Kovalenko, N. V., et al. (2020). A phase III, randomized, double-blind, multicenter study to compare the efficacy, safety, pharmacokinetics, and immunogenicity between SB8 (proposed bevacizumab biosimilar) and reference bevacizumab in patients with metastatic or recurrent nonsquamous non-small cell lung cancer. Lung cancer 146, 12–18. doi:10.1016/j.lungcan.2020.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Roche Registration GmbH (2005). Avastin European public assessment report. Available at: https://www.ema.europa.eu/en/documents/scientific-discussion/avastin-epar-scientific-discussion_en.pdf (Accessed August 12, 2023).

Google Scholar

Roche Registration GmbH (2023). Avastin summary of product characteristics. Available at: https://www.ema.europa.eu/en/documents/product-information/avastin-epar-product-information_en.pdf (Accessed August 12, 2023).

Google Scholar

Samsung Bioepis Co. Ltd (2023). Aybintio (bevacizumab for injection) product monograph. Available at: https://pdf.hres.ca/dpd_pm/00069363.PDF (Accessed August 12, 2023).

Google Scholar

Samsung Bioepis NL B.V. (2023a). Onbevzi summary of product characteristics. Available at: https://www.ema.europa.eu/en/documents/product-information/onbevzi-epar-product-information_en.pdf (Accessed August 12, 2023).

Google Scholar

Samsung Bioepis NL B.V. (2023b). Aybintio summary of product characteristics. Available at: https://www.ema.europa.eu/en/documents/product-information/aybintio-epar-product-information_en.pdf (Accessed August 29, 2023).

Google Scholar

Shepshelovich, D., Xu, W., Lu, L., Fares, A., Yang, P., Christiani, D., et al. (2019). Body mass index (BMI), BMI change, and overall survival in patients with sclc and NSCLC: a pooled analysis of the international lung cancer consortium. J. Thorac. Oncol. 14 (9), 1594–1607. doi:10.1016/j.jtho.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Shin, D., Lee, Y. J., Choi, J., Lee, D., Park, M., and Petkova, M. (2020). A phase I, randomized, single-dose pharmacokinetic study comparing sb8 (bevacizumab biosimilar) with reference bevacizumab in healthy volunteers. Cancer Chemother. Pharmacol. 86 (4), 567–575. doi:10.1007/s00280-020-04144-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Symeonides, S., and Gourley, C. (2015). Ovarian cancer molecular stratification and tumor heterogeneity: a necessity and a challenge. Front. Oncol. 5 (229), 229. doi:10.3389/fonc.2015.00229

PubMed Abstract | CrossRef Full Text | Google Scholar

Taiwan Food and Drug Administration (2018). Onbevzi license details. Available at: https://info.fda.gov.tw/MLMS/H0001D.aspx?Type=Lic&LicId=60001193 (Accessed August 14, 2023).

Google Scholar

Tang, S., Zhang, Y., Liang, W., and He, J. (2016). Predictive value of vascular endothelial growth factor A (VEGF-A) for bevacizumab-based treatments across advanced cancers: a meta-analysis based on eight phase III randomized control trials involving 4,523 patients. Oncol 27 (6), vi31. doi:10.1093/annonc/mdw363.49

CrossRef Full Text | Google Scholar

Wang, Y., Sung, C., Dartois, C., Ramchandani, R., Booth, B. P., Rock, E., et al. (2009). Elucidation of relationship between tumor size and survival in non-small cell lung cancer patients can aid early decision making in clinical drug development. Clin. Pharmacol. Ther. 86 (2), 167–174. doi:10.1038/clpt.2009.64

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y. C., Wang, Y., Schrieber, S. J., Earp, J., Thway, T. M., Huang, S. M., et al. (2019). Role of modeling and simulation in the development of novel and biosimilar therapeutic proteins. J. Pharm. Sci. 108 (1), 73–77. doi:10.1016/j.xphs.2018.10.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Weise, M., Bielsky, M. C., De Smet, K., Ehmann, F., Ekman, N., Giezen, T. J., et al. (2012). Biosimilars: what clinicians should know. Blood 120 (26), 5111–5117. doi:10.1182/blood-2012-04-425744

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., Zhang, S., Xu, T., Zhan, M., Chen, C., and Zhang, C. (2022). Efficacy and safety of bevacizumab biosimilars compared with reference biologics in advanced non-small cell lung cancer or metastatic colorectal cancer patients: a network meta-analysis. Front. Pharmacol. 13, 880090. doi:10.3389/fphar.2022.880090

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, P., Sy, S. K. B., and Skerjanec, A. (2018). Application of pharmacometric analysis in the design of clinical Pharmacology studies for biosimilar development. AAPS J. 20, 40. doi:10.1208/s12248-018-0196-7

PubMed Abstract | CrossRef Full Text | Google Scholar