Optimal Scheduling of Bevacizumab and Pemetrexed/Cisplatin Dosing in Non-Small Cell Lung Cancer.

Bevacizumab-pemetrexed/cisplatin (BEV-PEM/CIS) is a first-line therapeutic for advanced nonsquamous non-small cell lung cancer. Bevacizumab potentiates PEM/CIS cytotoxicity by inducing transient tumor vasculature normalization. BEV-PEM/CIS has a narrow therapeutic window. Therefore, it is an attractive target for administration schedule optimization. The present study leverages our previous work on BEV-PEM/CIS pharmacodynamic modeling in non-small cell lung cancer-bearing mice to estimate the optimal gap in the scheduling of sequential BEV-PEM/CIS. We predicted the optimal gap in BEV-PEM/CIS dosing to be 2.0 days in mice and 1.2 days in humans. Our simulations suggest that the efficacy loss in scheduling BEV-PEM/CIS at too great of a gap is much less than the efficacy loss in scheduling BEV-PEM/CIS at too short of a gap.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑Bevacizumab is currently recommended for concomitant administration with antiproliferative drugs. Several studies have indicated that sequential scheduling of bevacizumab and antiproliferatives is superior to concomitant administration. Precisely determining the optimal length of time between bevacizumab and antiproliferative administration has remained a challenging goal.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑What is the optimal schedule for the administration of bevacizumab and pemetrexed‐cisplatin in mice and in humans?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑This study more precisely estimates the optimal schedule for bevacizumab‐pemetrexed/cisplatin than previous studies. This study also scales the mathematical model used to make those predictions in humans instead of only mice.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

☑This study further develops a semimechanistic model that can be used to describe the effect of administering an antiangiogenic and antiproliferative sequentially. The parameter and interindividual variability estimates can be used in future related studies to improve future drug development. This study also demonstrates how this model can be scaled to make predictions in humans.

Bevacizumab‐pemetrexed/cisplatin (BEV‐PEM/CIS) combination therapy has been shown to be an effective first‐line and maintenance therapy for non‐small cell lung cancer (NSCLC) in phase II and phase III clinical trials.1, 2 PEM inhibits the enzymes necessary for pyrimidine and purine synthesis—primarily thymidylate synthase, which is necessary for thymidine synthesis and tumor cell replication.3 CIS is an alkylating agent that crosslinks adjacent N7 centers on purine residues, damaging DNA, disrupting repair, and disrupting purine synthesis.4, 5, 6 Disrupting DNA substrate supply results in S‐phase arrest, DNA repair disruption, and eventually apoptosis.7, 8 CIS also significantly disrupts calcium and reactive oxygen species regulation, inducing cellular lesions that further sensitizes cancer cells to apoptosis.6

In contrast to the effect of PEM/CIS, i.e., DNA damage, BEV is an anti–vascular endothelial growth factor humanized monoclonal antibody. Vascular endothelial growth factor is an angiogenic potentiator that promotes the growth of endothelial tissue necessary for arteries, veins, and lymphatics. By limiting neovascular growth, and therefore blood delivery to neoplasms, BEV exhibits limited antiproliferative properties.9

More importantly, BEV transiently induces a pruning effect on neovascular beds, which normalizes blood supply to neovascularly dense tissues (i.e., tumors).10, 11, 12 By normalizing blood supply, BEV enhances chemotherapeutic (i.e., PEM/CIS) delivery to neoplasms.13, 14

The effects of BEV‐PEM/CIS are generalized, i.e., any cell capable of uptaking the drugs are susceptible to their effects, especially rapidly dividing cells such as myeloid cells.15 Accordingly, BEV‐PEM/CIS has a narrow therapeutic window and generalized side effects.16 Previous studies on BEV‐PEM/CIS suggest that the sequential administration of BEV‐PEM/CIS (i.e., BEV before PEM/CIS) outperforms concomitant scheduling of BEV‐PEM/CIS in treating NSCLC.11, 17, 18, 19 This makes BEV‐PEM/CIS an attractive target for scheduling optimization via modeling and simulation, as a range of practical predictions—such as optimal scheduling in humans—can be made without the considerable time and resource investment required to conduct in vivo experiments. These predictions can be used to guide future studies, greatly accelerating drug development.20

In our previous work on BEV‐PEM/CIS published in Imbs et al.,17 mice with NSCLC tumor xenografts were administered BEV with PEM/CIS combination therapy in either concomitant or delayed (i.e., BEV before PEM/CIS) scheduling. The NSCLC tumors had been modified such that tumor growth could be tracked over time via either bioluminescence or fluorescence. Following previous theoretical investigations, the data set generated from the mice with bioluminescent tumors was used to develop a semimechanistic pharmacokinetic (PK)/pharmacodynamic (PD) model for tumor dynamics in response to BEV‐PEM/CIS.21, 22 The model was then used to predict the optimal scheduling gap between BEV and PEM/CIS administration.

The aim of this follow‐up modeling work was to both refine and expand on previous results on BEV‐PEM/CIS combination therapy using the much larger fluorescence data set generated in Imbs et al.17 We first showed that the semimechanistic model previously developed better explained the data than comparable models (i.e., we validated the previously developed structural model). We then refined the parameter estimates of the model and used it to predict the optimal scheduling gap between BEV and PEM/CIS administration. Next, we used stochastic simulations to explore the marginal loss in therapeutic efficacy when BEV‐PEM/CIS was administered at a suboptimal gap and the effect of BEV dose scaling on population optimal gap as well as the interindividual variability (IIV) of the optimal gap. Here, marginal loss in therapeutic efficacy refers to the loss in treatment efficacy from either scheduling the drugs with a gap shorter or longer than the optimal gap, e.g., 1 day earlier or later than the optimal gap. Lastly, using human PK/PD models and parameter estimates from the literature, we were able to scale the model to estimate the optimal scheduling of BEV‐PEM/CIS in humans.

Methods

Experimental procedure

Comprehensive details on the animals and the experimental procedure are available in Imbs et al.17 Briefly, on day 0 of the experiment, tumors (ca. 120,000 cells) consisting of H460 human NSCLC transfected with luciferase and the tdTomato gene (H460 Luc+ tdTomato+; Perkin Elmer France) were injected ectopically into the left flank of 90 mice. The animals were pathogen‐free, immunocompromised, 6‐week‐old, female Swiss nude mice (Charles River, France). The mice were randomized into one of five treatment groups. The first study group (control) received no treatment. The second treatment group (PEM/CIS) was administered both 100 mg/kg of PEM intraperitoneally (i.p.) and 3 mg/kg of CIS i.p. on days 14, 28, and 42 of the experiment. The third, fourth, and fifth treatment groups (BEV‐PEM/CIS) received the same PEM/CIS treatment as the second experimental group.

In addition, the BEV‐PEM/CIS treatment groups were administered 20 mg/kg i.p. of BEV either concomitantly with the PEM/CIS administrations (group 3), 3 days prior to each PEM/CIS administration (group 4), or 8 days prior to each PEM/CIS administration (group 5); see Table for administration tabulation.

Tumor growth was monitored on a minimum biweekly basis using Ivis Spectrum imager (Perkin Elmer France), and the images were acquired and analyzed using the Living Image 6.0. software (Perkin Elmer France).

The mice were provided paracetamol‐supplemented water (e.g., 80 mg/kg/day) to prevent disease‐related pain. Animals showing signs of distress, pain, cachexia (i.e., loss of 10% of body mass), or extensive tumor proliferation (i.e., within 2–3 cm) were euthanized. All animals were euthanized on day 87 of the experiment. All experiments were approved by the local ethical committee at French Ministère de l'Education Nationale, de l'Enseignement Supérieur et de la Recherche, and registered as 2015110616255292.

PK/PD structural model building and evaluation

The PK models for BEV, PEM, and CIS were derived from previously published PK models in mice.23, 24, 25 The parameters for these models were fixed to the typical values from those studies and assumed no IIV.

The PD model was selected from a series of sequentially fit tumor growth and drug effect models. First, using only the control data set, the exponential, linear‐exponential, and Gompertz growth models were cross‐evaluated as models of unperturbed tumor growth.26 Then, incorporating the full data set into the fit, the log‐kill effects of BEV, PEM, and CIS were each considered. The interaction effect between BEV and PEM/CIS was included to represent the synergistic effect of BEV. Following previous work for the effect of cytotoxic drugs, three cellular death compartments were included in the PD portion of the model to represent the delay between cellular damage as a result of PEM/CIS and cell death.27

Competing models were evaluated numerically using Bayesian information criteria and the precision of parameter estimates—defined as the relative standard error of the estimate (RSE). Observed vs. predicted plots, individual fit plots, and visual predictive checks (VPCs) were produced to graphically assist model evaluation (as automated in Monolix 2018R2). VPCs were produced using the default estimation process for VPCs as of Monolix 2018R2, i.e., to create the 90% prediction intervals for the 10th, 50th, and 90th percentiles, 500 simulations are performed using random individual parameters and the design structure of the experiment.

After model selection, the statistical correlations between random effects were explored via visual inspection. Correlations plots between random effects were produced, defined as ηi,t,ϕ1 vs ηi,t,ϕ2 i.e. the random effect η, of individual i, at time t, of parameter 1, i.e. ϕ1, vs the random effect η, of individual i, at time t, of parameter 2, i.e. ϕ2. The full posterior distribution of the parameters were used in place of empirical Bayes estimates to avoid visual artifacts as a result of shrinkage as suggested by Lavielle and Ribba28 and Pelligand et al.29 Statistical correlations between random effects were also numerically assessed using a Pearson correlation test at a P < 0.05 threshold.

Stochastic approximation expectation method (SAEM) convergence and final model parameterization were graphically assessed by an inspection of search stability, distribution of the individual parameters, distribution of the random effects, individual prediction vs. observation, individual fits, and distributions of the weighted residuals as well as VPCs.

The precision of parameter estimates was numerically assessed using RSE. The normality of random effects distributions, the normality of individual parameter distributions, and the normality of the distribution of residuals were each numerically assessed using a Shapiro–Wilk test (P < 0.05). The centering of the distribution of residuals (i.e., centered on zero) was numerically assessed using a Van Der Waerden test (P < 0.05).

Parameter stability was assessed by comparing parameterizations resulting from random initial starting value selection—as implemented in the Monolix assessment suite. The assessment suite performs 5 SAEM parameterizations in series using random initial parameter values uniformly drawn from the interval from approximately 60–160% of final parameter estimates. The SAEM of the individual parameterizations was then tracked between runs, giving a range of parameter value estimates, RSEs, and log‐likelihoods to compare. This assessment was used to ensure that the algorithm did not converge to a log‐likelihood local minimum during the process of producing final parameter estimates. The settings described are the default as of Monolix 2018R2.

Simulations

Simulations were performed in R 3.4.4 using Simulx 3.3.030 to simulate from Monolix run files. First, a function was built that accepted treatment schedule, parameter substitutions, dose, and number of individuals as input and produced a simulated population as an output. This function was simply a convenience wrapper of Simulx for automation purposes and was verified by reproducing VPCs per treatment group. Simulation Set 1 was used to predict the optimal gap between administration of BEV and PEM/CIS. Simulation Set 2 produced an estimate of the IIV of the optimal gap. Simulation Set 3 examined the anticipated effect of varying the dose of BEV on the optimal gap. Simulation Set 4 scaled predictions of BEV‐PEM/CIS efficacy to humans. All simulations were population‐level responses (i.e., simulated without RSE or IIV) except for Simulation Set 2 (simulated without RSE and with IIV), which used 1,000 Monte Carlo samples. Further details are provided in the .

Quality assessment

All mlxtran and R codes were assessed for quality control by an independent evaluator.

Results

Error model

Measurement error was best described using a log‐normal, constant‐error model (Eq.  (1)). The natural log of each individual measurement, ln(y ij) with individual i and repetition j, was modeled as a measurement centered on the natural‐log mean of y ij over j, i.e., ln, in addition to some residual error, ϵij , normally distributed, centered on zero, and with standard deviation a.

The standard goodness‐of‐fit graphics and numerical analyses supported that a log‐constant error model best fit the data. The log‐constant error model was not rejected for both the Kolmogorov–Smirnov test for normality (P = 0.22) nor the Pearson Χ2 normality test (P = 0.0509). In contrast, a constant error model was rejected by these two tests (see Figure for further details).

PK/PD structural model building

No outliers were identified during the initial data exploration. Therefore, no collected data were excluded from model building.

The PK of BEV, PEM, and CIS were each modeled using one‐compartment models with first‐order i.p. absorption and first‐order elimination based on literature descriptions and PK parameter estimates (Table ).23, 24, 25 Random effects (i.e., ηpk) were set to 0 as individual PK was not reported in these experiments.

A Gompertz function (Eq.  (2)) was found to best describe the unperturbed tumor growth V(t) based on its fit performance over competing models, low RSE on parameter estimates, and literature‐established descriptive quality.

Because of the relative sparseness in sampling, fitting more complex semi mechanistic models of growth incorporating cellular quiescence, biphasic growth, and cell‐type heterogeneity did not produce model fits with lower Bayesian information criteria or comparable precision in parameter estimates (RSE) to the simpler Gompertz model of tumor growth. Parameters α and β represent the proliferation rate of the tumor cells and the rate of exponential decrease of the tumor relative growth rate, respectively. V c is the unit value of relative fluorescence units (RFU) corresponding to one cell, i.e., the proportionality constant between RFU and the number of cancer cells in the fluorescent volume.

V c was estimated externally from Monolix by conducting a naive‐pooled, linear regression on the natural log of the full data set. The regression gave a rough estimate of V 0, which was then scaled by the approximate number of cells injected at time 0 (ca. 120,000 cells) to derive V c = 5.064 × 10−4 RFU.

After selecting an appropriate growth model, the log‐kill effects of PEM, CIS, and BEV were each considered in parallel. The log‐kill effect of BEV was estimated as insignificant and removed from the model. The estimation of the log‐kill effects of PEM and CIS were found to be highly correlated. To reduce model complexity, only their combined concentration, C(t), and a corresponding log‐kill parameter, γ, were considered in the final model (Eq.  (3)).

A(t) represents the plasma concentration of BEV. Q(t) represents the synergistic effect of improved vascular quality. In brief, the increase in neoplasm vascular quality because of BEV typically occurs within a period of a few days after administration. To represent this delay in effect, time (t) was delayed by τ. Parameter δ represents the proportional increase in PEM/CIS efficacy as result of vascular quality improvement under BEV therapy.

The estimation of τ was bounded between 0 and 10 using the link function τ = τbound·10, where logit(τbound) ~ N(0,1). All other parameters were best estimated as log‐normally distributed. The full statistical representation of individual parameters, ϕi, estimated via SAEM is shown in Eq.  (4), where the full structural model is denoted by F.

Cellular death as a result of chemotherapeutic treatment was modeled as a three‐compartment transition from the growth compartment to death.27 The compartments are labeled Z 1, Z 2, and Z 3 (numbering respective to their order), and the transition between compartments is governed by the intercompartmental clearance parameter k.

After a period of manual exploration, k was set to the value of 0.3. This choice is consistent with the parameterization made in Imbs et al.17 This choice also limits the total transition time from the tumor mass compartment to cellular death to the order of a day, which is consistent with the upper limits of cellular death clearance.31

The full tumor size, N, was the sum of the size of unperturbed cells, V, as well as the size of damaged cells undergoing cellular death, i.e., Z 1 + Z 2 + Z 3.

No correlations between random effects were statistically significant enough to be included in the final model (Figure ). Full parameter estimates and the model diagram are provided in Table  and Figure  respectively. Model diagnostics are collected in Figures  , , .

Table 1

Pharmacodynamic model parameters for tumor proliferation in non‐small cell lung cancer–xenografted mice

Model parameter	Estimate	SE	RSE (%)	IIV	IIV (CV%)
α	0.77 day−1	0.0081	1.06	0.037	4.76
β	0.04 day−1	0.0005	1.16	0.015	33.55
γ	35.74 (mg day)−1	6.36	17.79	0.46	1.28
δ	3.73 (mg day)−1	0.92	24.70	0.35	9.51
τ	1.19 days	0.017	1.50	0.17	14.50
V 0	3.4 RFU	0.27	22.55	1.73	143.03
k	0.30 day−1	—	—	—	—
Residual error variance	0.43	0.01	2.55	—	—

Model parameter estimates (fixed and random effects) as well as standard errors as determined by the stochastic approximation expectation‐maximization algorithm as implemented in Monolix 2018R2. α and β represent the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate, respectively. γ was used to model the log‐kill effect of pemetrexed and cisplatin's combined concentration. δ represents the proportional increase in pemetrexed/cisplatin efficacy as a result of vascular quality improvement under bevacizumab therapy. τ was the time delay in bevacizumab effect. V 0 represents the volume of cells at time 0. k governs intercompartmental clearance between cellular death compartments. CV, coefficient of variation; IIV, interindividual variability; RFU, relative fluorescence unit; RSE, relative standard error of the estimate.

Figure 1

Structural model diagram. The scheme of the structural model is depicted to the right. Unperturbed cells grow at rate governed by α and β. When a cytotoxic is introduced into the system, the cytotoxic impairs the growth of the tumor by sending cells into a death succession. The parameter that determines the cytotoxic efficacy, γ, is scaled by both the concentration of cytotoxics, C(t), and the volume of the tumor, V(t). Bevacizumab improves vascular quality, Q(t), after time delay, τ, which scales the cytotoxic effect by parameter δ. When a cell is damaged by cytotoxics, it begins a progression from unperturbed growth—compartment V(t)—to damage compartments Z 1 through Z 1. Eventually the cell exits the tumor volume as it dies. The rate of transfer between damage compartments is governed by intercompartmental clearance parameter k.

Figure 2

Standard goodness‐of‐fit diagnostic plots. On the left is individual predictions vs. observations and on the right are the individualized weighted residuals (IWRES) vs. time. During model fitting, observations were natural log–transformed to stabilize predictions. Therefore, residuals, predictions, and observations are natural log–transformed in these figures. The predictions are approximately normally distributed. On the left, the one‐to‐one prediction line is the center solid black line, the spline (average agreement between individual prediction and observation) is solid orange, and the dashed black lines are the borders of the 90% prediction interval. On the right, the zero residual error line is the center dashed black line and the spline is solid orange. The dashed black lines are the borders of the 90% prediction interval.

Figure 3

Sample of individual fits. The blue dots represent individual observations, whereas the solid violet line represents individual fits. Some mice—e.g., mice in the control group—were not tracked for the full course of the study. These animals experienced complications as a result of the experiment and either spontaneously passed or were euthanized to prevent excessive suffering. PEM‐CIS, pemetrexed‐cisplatin; RFU, relative fluorescence unit.

Figure 4

Combined and stratified visual predictive checks. The blue lines are the 10th, 50th, and 90th empirical percentiles calculated for each unique value of time. Blue and pink areas represent 90% prediction intervals for the 10th (blue), 50th (pink), and 90th (blue) percentiles. Prediction intervals are calculated by Monte Carlo simulation. To create prediction intervals for each unique value of time, 500 simulations are performed using random individual parameters. The red areas and red‐circled points represent areas where empirical measurements fall outside of the bounds of the 90% prediction intervals. PEM‐CIS, pemetrexed‐cisplatin; RFU, relative fluorescence unit; VPC, visual predictive check.

Mouse simulations

Simulating the experimental treatments with a range of administration gaps from 0–10 days (step size = 0.1 day) suggested that the optimal time delay between scheduling BEV and PEM/CIS in mice is 2.0 days (Simulation Set 1).

The simulated IIV of the optimal gap was relatively small. Only three values of individual optimal gap were produced. Of the virtual animals, 96.5% had an individual optimal gap of 2.0 days, 1.0% of the virtual animals had an individual optimal gap of 2.1 days, and 2.5% of the virtual animals had an individual gap of 1.9 days (Simulation Set 2).

Scaling the dosage of BEV to either 30 mg/kg or 10 mg/kg produced no effect in the estimated optimal gap and produced no effect in the IIV of the optimal gap (Simulation Set 3).

Human simulations

Simulations of the typical human response to chemotherapy and BEV were performed using IV administration, two‐compartment absorption, and first‐order elimination models and parameters.32, 33, 34 Dosage, infusion time, and frequency of administration for BEV‐PEM/CIS were adapted from the induction phase of the AVAPERL phase III clinical trial in BEV‐PEM/CIS.2 The average adult weight and body surface area were obtained from the Centers for Disease Control and Prevention and literature estimates, respectively.35, 36

Except for the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate (i.e., α and β), the PD model and parameterization (γ, δ, τ, κ) were reused exactly as they were determined in the mouse portion of the model. The α and β estimates were obtained from Bilous et al.,37 where clinical NSCLC doubling times reported in Friberg and Mattson38 were used to estimate population α and β for NSCLC in humans. The value of V c came from the classical assumption that a 1‐mm3 volume of tumor cells is approximately 106 cells.39 V0 was arbitrarily set to 3 cm3. This scaling procedure resulted in an adapted model primarily constructed to predict optimal gap in humans. In contrast, the model relies on several implicit assumptions for scaling tumor response (see ).

The full PK/PD model was then used to simulate the typical cancer growth under various administration schedules with a starting tumor volume of 3 cm3. Parameter estimates are reported in Table  , the simulation administration schedule is reported in Table S2, and the simulation summaries are depicted in Figure  . The estimated optimal gap between BEV and PEM/CIS administration in humans was 1.2 days.

Table 2

Model parameterization for simulations of non‐small cell lung cancer treated with bevacizumab‐pemetrexed/cisplatin in humans

Drug	Pharmacokinetics
	Pemetrexed	Cisplatin	Bevacizumab
V 1 (L)	12.9	22.3	2.8
V 2 (L)	3.38	77.0	2.9
Q (L/day)	20.7	456.0	0.6
CL (L/day)	131.9	6.5	0.2
Infusion time (minute)	10.0	120.0	60.0

Except for the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate (i.e., α and β), the pharmacodynamic model and parameterization were reused exactly as they were determined in the mouse portion of the model. α and β estimates were obtained from Bilous et al., 37 where clinical non‐small cell lung cancer doubling times reported in Friberg and Mattson38 were used to estimate population α and β for non‐small cell lung cancer in humans. The value of V c came from the classical assumption that a 1‐mm3 volume of tumor cells is approximately 106 cells.39 V 0 was arbitrarily set to 3 cm3. V 1, V 2, Q, and CL represent volume of compartment one and peripheral compartment two, intercompartmental clearance, and clearance from compartment one, respectively.

Figure 5

Human pharmacodynamic and pharmacodynamic simulations summary. To produce these figures, bevacizumab was administered anywhere from 0–10 days (in steps of 0.1) before pemetrexed/cisplatin was administered. Tumor growth was simulated from 0–67 days with no interindividual variability and no relative standard error of the estimate. In the top figure, the AUC of tumor growth vs. gap (0–10 days) is depicted. In the middle figure, tumor dynamics over time, with gap indicated by color, are depicted. In the bottom panel, the pharmacokinetics of bevacizumab–pemetrexed/cisplatin is depicted with gap indicated by color. The top figure indicates that the optimal scheduling gap is 1.2 days and the middle figure depicts the difference in tumor volume between administration gaps. The patient with optimal scheduling had a final tumor volume ~30% the size of the concomitant scheduling.

Discussion

By normalizing tumor vasculature, BEV improves the delivery of PEM/CIS to tumors and consequently PEM/CIS efficacy. PEM and CIS each have a narrow therapeutic window and high toxicity. It is therefore critical that BEV‐PEM/CIS doses are administered as efficiently as possible. This makes BEV‐PEM/CIS a natural fit for modeling and simulation studies, as the drug scheduling can be optimized without the need for multiple time and resource intensive in vivo studies. In this analysis, we conducted an in silico study of the optimal administration of BEV‐PEM/CIS in a xenograft and human model of NSCLC by constructing a mathematical model of tumor dynamics in response to BEV‐PEM/CIS. In constructing that model, we were able to validate and refine previous modeling in BEV‐PEM/CIS. Greater precision in parameter estimates was achieved through the external estimation of V c and the external validation of the residual error model as well as by using the larger fluorescence data set to obtain final parameter estimates. Then, after exploring a range of predictions in mice, we scaled our model to predict optimal scheduling in humans.

In the error‐modeling portion of the experiment, we demonstrated a strategy through which the choice of the error model can be validated externally to the primary data set by including supplementary data collection in the experimental design. This simplified the error‐modeling step in the model‐building process.

The next stage of this study consisted of determining whether the semimechanistic model of tumor dynamics in response to BEV‐PEM/CIS developed in Imbs et al.17 best fit the unfit fluorescence data from the same study. During model building, we attempted to balance our model‐building procedure between model performance (empirical fit) and the underlying biology, an approach often referred to as the middle‐out approach.40, 41

In selecting potential PD models of tumor growth, several semimechanistic models were fit to the experimental data. The Gompertz model and linear exponential model performed comparably. The parameters of the Gompertz model were estimated with greater precision than the parameters of the linear exponential model (RSE), where the linear exponential model was fit with lower Bayesian information criteria than the Gompertz model. Ultimately, the Gompertz model was chosen over the linear exponential model because of the physiological relevance of its construction.

The parameterization of the final model was slightly unstable because of modest overparameterization. To compensate for this, k was fixed to a reasonable physiological estimate to improve the precision of parameter estimates, and the search for τ was bounded to reduce spurious individual parameter estimates.

During structural model building, we confirmed the validity of the model previously published in Imbs et al.17 We also reconfirmed the efficacy improvement of BEV‐PEM/CIS dosing over PEM/CIS or control. We observed that a 3‐day gap in scheduling is superior to both concomitant scheduling and an 8‐day gap in scheduling. We were also able to build on previous work by identifying with greater precision the parameters underlying the mathematical model of BEV‐PEM/CIS in NSCLC–tumor bearing mice.

In our mouse simulations, the final tumor volume (after 67 days) in the optimal scheduling group with BEV‐PEM/CIS (gap = 2.0 days) was 88.5% of the size of final tumor volume in the concomitant scheduling group. This is consistent with our experimental results, i.e., that mice administered BEV approximately 2.0 days before PEM/CIS have a moderately better response (i.e., greater tumor size reduction) to BEV‐PEM/CIS than mice who are administered BEV‐PEM/CIS concomitantly. We also found, through simulation, that scaling the dose of BEV had no effect on the optimal gap and that IIV on gap is low.

Predictions made by our model agree with previous findings in BEV‐PEM/CIS scheduling. The order of the optimal scheduling delay (2.0 days) is within the 1–5 day gap predictions of previous studies.11, 18, 19 Studies in tumor perfusion and BEV showed day 1 and day 4 decreases in tumor perfusion, which is consistent with the marginal predictions in our study, i.e., optimal perfusion should be on the order of 2 days with comparable marginal losses on either side of that minimum.42

After adapting our model to make simulations in humans, we predicted a robust response to both sequential and concomitant BEV‐PEM/CIS scheduling. The final tumor volume (after 85 days) in the optimal scheduling group (gap = 1.2 days) was 52% of the size of the final tumor volume in the concomitant scheduling group. If these predictions are accurate, scheduling optimization could result in significant improvement in BEV‐PEM/CIS combination therapy efficacy with no increase in toxicity.

Because the dose of the therapeutics, their PK, the administration method, and dosing schedules are different in humans vs. tumor‐bearing mice, a shift in the optimal gap is expected between species. In this study, the discrepancy between the optimal gap in humans vs. NSCLC tumor‐bearing mice is driven by the slow i.p. absorption of BEV in mice (Table ).

When exploring marginal efficacy loss in suboptimal administration schedules, we consistently found that the marginal cost of scheduling BEV and PEM/CIS too close together in time was greater than the marginal cost of scheduling BEV and PEM/CIS with too great of a gap in administration—in both mice and humans. This indicates that any potential clinical studies in antiangiogenics and cytotoxics should weight scheduling recommendations toward scheduling at slightly too large of gap.

Lastly, we would like to acknowledge some limitations of this study. The tumor microenvironment is known to be complex and varied. Tumor tissues contain necrotic pockets, heterogenous and dynamic microvasculature, and various submutations that result in differential local growth rates and drug sensitivities. Furthermore, BEV efficacy is potentially disease‐state dependent. Considering this biological heterogeneity and identifying biomarkers for measuring individual response would greatly improve future model predictions and model scalability between species as well as modeling‐based individualized therapy development.

In summary, our analysis confirms previous findings in BEV‐PEM/CIS scheduling while improving precision of parameter estimates, improving prediction quality and detail, and scaling the model to predict the optimal scheduling of BEV‐PEM/CIS in humans. Antiangiogenics will continue to be useful agents in oncology. There are currently several other antiangiogenics regularly used in combination with cytotoxics that could potentially benefit from sequential administration (i.e., antiangiogenic then cytotoxic).43 Of note, BEV is currently only approved for concomitant administration with chemotherapy in all of its indications, e.g., lung cancer, breast cancer, gastric cancer, and so on. This contrasts with the optimized sequential scheduling that model simulations suggest.

There is a recent trend to develop model‐informed drug development to optimize anticancer therapy. Our work highlights how mathematical modeling could help to refine clinical treatment modalities. The semimechanistic nature of this model allows it to be modularly reconfigured to extend predictions to other antiangiogenics as well as novel therapeutic paradigms such as the immune checkpoint inhibitor monoclonal antibody pembrolizumab.44 This work continues to lay the foundation for building systems pharmacology models for the effect of antiangiogenic and antiproliferative combination therapy in advanced NSCLC. Tortuous vasculature is a phenotype exhibited by many solid tumors, and predicting optimal antiangiogenic scheduling could greatly increase the efficacy of future oncology therapeutics and combination therapies.45

Funding

No funding was received for this work.

Conflict of Interest

J.C. received fees for participating in the board meetings for Roche Laboratory. All of the other authors declared no competing interests for this work.

Author Contributions

B.K.S., S.B., and J.M. wrote the manuscript. J.C., F.B., and S.B. designed the research. A.B. and J.C. performed the research. B.K.S., S.B., and J.M. analyzed the data.

Supplementary Figures.

Click here for additional data file.

Supplementary Methods.

Click here for additional data file.

Project Code—PKPD Model of Beva–Cyto R01.

Click here for additional data file.