# Optimizing First-Line Therapeutics in Non-Small Cell Lung Cancer: Insights From Joint Modeling and Large-Scale Data Analysis

Benjamin K Schneider, Sébastien Benzekry, Jonathan P Mochel

► To cite this version:

Benjamin K Schneider, Sébastien Benzekry, Jonathan P Mochel. Optimizing First-Line Therapeutics in Non-Small Cell Lung Cancer: Insights From Joint Modeling and Large-Scale Data Analysis. CPT: Pharmacometrics and Systems Pharmacology, 2025, 10.1002/psp4.70079. hal-05198474

HAL Id: hal-05198474

https://inria.hal.science/hal-05198474v1

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE OPEN ACCESS

# Optimizing First-Line Therapeutics in Non-Small Cell Lung Cancer: Insights From Joint Modeling and Large-Scale Data Analysis

Benjamin K. Schneider¹ ID | Sebastien Benzekry¹,² ID | Jonathan P. Mochel¹,³ ID

¹Iowa State University College of Veterinary Medicine, Ames, Iowa, USA | ²Computational Pharmacology and Clinical Oncology Department, Inria Sophia Antipolis - Méditerranée, Cancer Research Center of Marseille, Inserm UMR1068, CNRS UMR7258, Aix Marseille University UM105, Marseille, France | ³Precision One Health Initiative, Systems Pharmacology, University of Georgia, Athens, Georgia, USA

**Correspondence:** Benjamin K. Schneider (schnebeni@gmail.com)

**Received:** 10 November 2024 | **Revised:** 10 June 2025 | **Accepted:** 23 June 2025

**Funding:** Research was funded primarily by Iowa State University, as well as a salary supplement via Ceva Sante Animale.

**Keywords:** bevacizumab | lung cancer | mathematical modeling | methodology | model based meta analysis | non-linear mixed effects | non-small cell lung cancer | NSCLC | population pharmacokinetics-pharmacodynamics | resistance

## ABSTRACT

Non-small cell lung cancer (NSCLC) is often intrinsically resistant to several first- and second-line therapeutics and can rapidly acquire further resistance after a patient begins treatment. Treatment outcomes are, therefore, significantly impacted by the optimization of scheduling. Previous preclinical research has suggested that scheduling bevacizumab sequentially with combination antiproliferatives could improve clinical outcomes. Mathematical modeling is a well-suited tool for investigating this proposed modification. To address this critical need, individual patient tumor data from 11 clinical trials in NSCLC have been collated and used to develop a semi-mechanistic model of NSCLC growth and response to the therapeutics represented in those trials. Precise estimates of clinical parameters fundamental to cancer modeling have been produced, such as the rate of acquired resistance to various pharmaceuticals, the relationship between drug concentration and cancer cell death, as well as the fine dynamics of vascular remodeling in response to bevacizumab. In a reserved portion of the dataset, this model predicted the efficacy of individual treatment time courses with an average difference between final prediction and observation of 59.7% after a single tumor measurement and 11.7% after three successive tumor measurements. A delay of 9.6 h between pemetrexed-cisplatin and bevacizumab administration is predicted to optimize the benefit of sequential administration. At this gap, approximately 93.5% of simulated patients benefited from a gap in administration compared with concomitant administration. Of those simulated patients, the mean improvement in tumor reduction was 20.7%. This suggests that scheduling a modest gap between the administration of bevacizumab and partner antiproliferatives could meaningfully improve patient outcomes in NSCLC.

## 1 | Introduction

Lung cancer is the leading cause of cancer mortality nationwide, with an estimated 125,070 deaths expected in the United States

this last year 2024 [1]. Approximately 85% of those deaths are attributable to non-small cell lung cancer (NSCLC). The need for pharmacological intervention is dependent on the disease stage at the time of diagnosis. If NSCLC is detected early, surgical

Sebastien Benzekry and Jonathan P. Mochel co-last authors.

This is an open access article under the terms of the [Creative Commons Attribution-NonCommercial License](https://creativecommons.org/licenses/by-nc), which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

© 2025 The Author(s). CPT: Pharmacometrics & Systems Pharmacology published by Wiley Periodicals LLC on behalf of American Society for Clinical Pharmacology and Therapeutics.
## Summary

* What is the current knowledge on the topic?

* Non-small cell lung cancer (NSCLC) is a leading cause of mortality. Many treatments for NSCLC rely on combination therapy with bevacizumab, an anti-VEGF adjunct. Previous preclinical research has suggested that scheduling bevacizumab in sequence with combination antiproliferatives could significantly improve clinical outcomes.

* What question did this study address?

* To address this hypothesis, individual patient tumor data were collated from 11 clinical trials in NSCLC and used to develop a semi-mechanistic model of NSCLC growth and response to the various therapeutics, in combination with bevacizumab, represented in those trials.

* What does this study add to our knowledge?

* This study has produced precise estimates of several pharmacometric parameters fundamental to cancer modeling. Simulations from this model suggest that a delay of 9.6h between pemetrexed-cisplatin and bevacizumab administration optimizes the benefit of sequential administration. At this gap, approximately 93.5% of simulated patients benefited from a gap in sequential administration compared with concomitant administration. Of those simulated patients, the mean improvement in tumor reduction was 20.7%.

* How might this change drug discovery, development, and/or therapeutics?

* This review of clinical data suggests that scheduling a modest gap between the administration of bevacizumab and its partner antiproliferatives could meaningfully improve patient outcomes in NSCLC. This work indicates that future trials in NSCLC with bevacizumab should test staggered sequential scheduling. Furthermore, a model fit on a large human clinical trial dataset is provided for future modeling and simulation work in NSCLC.

resection followed by limited adjuvant therapy is generally successful in achieving complete remission [2, 3]. Late-stage NSCLC metastases are often too diffusely spread for surgery or radiotherapy to be effective; therefore, the occurrence of a diagnosis at an advanced stage necessitates treatment strategies predominantly centered on pharmacological interventions. The 5-year survival rate for patients diagnosed with localized NSCLC is approximately 58%. However, most patients with NSCLC are first diagnosed at the late stage of the disease, with a 5-year survival rate estimated at 6% in the presence of distant metastases [4, 5]. Ultimately, most patients with mid-to-late-stage NSCLC will be prescribed a platinum-based combination chemotherapy and/or immunotherapy, e.g., bevacizumab-carboplatin-paclitaxel, pembrolizumab-carboplatin-pemetrexed [6-8]. The most cutting-edge first-line combination therapies in NSCLC are atezolizumab-bevacizumab-carboplatin-paclitaxel combination therapy (established in IMpower150 study) [9] and pembrolizumab plus pemetrexed-platinum combination therapy (established in KEYNOTE-189) [10, 11]. Patients with advanced cases of NSCLC are often cycled through several antiproliferative

regimens due to high rates of acquired and intrinsic multidrug resistance [12–15].

Acquired or intrinsic drug resistance is a major cause of first-line therapeutic failure in NSCLC. In a previous study of in vitro, resected NSCLC by d'Amato et al., intermediate to extreme intrinsic resistance to carboplatin was found in 68% of NSCLC cultures (vs. 63% and 40% for cisplatin and paclitaxel, respectively) [12]. Likewise, in the recent KEYNOTE-001 study, NSCLC patients receiving pembrolizumab had an objective response rate of 19.4%, indicating that a vast majority of individuals did not significantly respond to therapy [15]. Finally, in the most recent ARIES observational cohort study of first-line NSCLC treatment involving bevacizumab, the failure rate for 1967 patients was approximately 51% [16]. Taken together, there is a *clear and critical unmet clinical need* to improve patient response to treatment and clinical outcomes in NSCLC.

Improving clinical outcome in NSCLC with existing therapeutics requires addressing at least three prominent challenges: (1). The disease is often late stage at the time of diagnosis, and treatment relies heavily on complex combination therapies. However, oncologists currently have few evidence-based tools to individualize and optimize scheduling of those therapies; (2). The success rate of first-line therapy is typically low, therapeutic dosing (amount, interval, gap) strategies are often unoptimized, and resistance to treatment is rapidly acquired [6–8]. (3). The window between the threshold for therapeutic efficacy and drug overdose can be narrow, especially with the relatively long elimination half-life of most antiproliferatives. Meeting these challenges is made even more difficult by the relatively sparse and sporadic collection of pharmacokinetic and pharmacodynamic data in clinical studies involving cancer.

The status quo as it pertains to therapeutic management and dosing decisions in NSCLC is that medications are administered according to recommended guidelines, and those guidelines are, in turn, established through one-to-one human clinical trials. A major limitation in this approach is that it mainly focuses on the *average* patient, such that clinicians have little evidence-based guidance to help them *individualize* treatment in patients who are unresponsive or achieve a partial response to therapy.

NLME models are a state-of-the-art tool for building mathematical descriptions of therapeutics from sparse data and continuously incorporating new information into predictions made via simulation. These simulations can be used to simulate experimental designs and guide treatment decisions, such as: *What happens if dosages are halved, but administrations are made twice as often? How should dosages of concurrent chemotherapeutics be modified? What is the ideal time after the first dose to start looking for signs of acquired resistance?* There are several paradigms for the development of NLME models. In this study, a semi-mechanistic modeling approach has been adopted to balance any limitations in the datasets with precise mathematical descriptions of NSCLC pathology developed in many previous studies. Semi-mechanistic models attempt to balance the desire for biologically relevant mechanisms with limitations in the datasets by modeling any prominent pathology with generalized empirical formulations [17, 18].
Before this study, a preliminary semi-mechanistic model of NSCLC growth dynamics was developed in a preclinical experimental setting [19]. Scaling this model to human biology, it was predicted that scheduling bevacizumab administration 1.2 days before pemetrexed-cisplatin administration has the potential to significantly improve tumor reduction (> 50%) as opposed to concomitant dosing [20]. In this follow-up study, the aim is to develop a predictive mathematical framework of first- and second-line therapeutic drugs for NSCLC by iterating on previous modeling efforts using a large (~11,000 patients) clinical trial dataset made available through the Vivli.org database (Table 1) [21].

The overall objectives in the present analysis were to build a model that (1) generalizes a previous model of NSCLC growth and response to bevacizumab-pemetrexed-cisplatin to the greater set of combination therapies and modes of action represented in this large clinical database, and (2) characterizes the time-course of resistance for those currently registered therapeutics for NSCLC.

In building this model, several empirical approximations of clinical parameters fundamental to cancer modeling have been estimated, such as the rate of acquired resistance to various pharmaceuticals, the relationship between drug concentration and cancer cell death, as well as the fine temporal dynamics of

anti-VEGF therapy. Then, using this high-powered model, an up-to-date prediction of optimal bevacizumab-pemetrexed-cisplatin was made in patients with NSCLC.

## 2 | Materials and Methods

See supplemental methods (Data S1) for sections on data processing and collation, non-linear mixed effects modeling, individual variation, and model evaluation.

### 2.1 | Literature Search and Review

To build this model, access to data was sought using very broad criteria. Data was prioritized from clinical trials where bevacizumab had been used in combination with other therapeutics to treat NSCLC. To build on previous modeling efforts, studies needed to include records of individual patient tumor sizes over time. After reviewing several potential data custodians to partner with, an application was sent for data access through the Clinical Study Data Request (CSDR) portal. Eleven different studies were identified that were available through CSDR's platform and met this study's requirements (Table 1) and access was permitted to a secure server containing the datasets beginning on January 6th, 2020.

<table><caption>TABLE 1 | Clinical data details.</caption><thead><tr><th>Trial registry ID</th><th>Phase</th><th>N</th><th>Therapeutics</th></tr></thead><tbody><tr><td>CSDR &amp; Vivli</td><td></td><td></td><td></td></tr><tr><td>NCT00800202</td><td>2</td><td>91</td><td>Bevacizumab, carboplatin, erlotinib, paclitaxel</td></tr><tr><td>AVF0757G</td><td>2</td><td>99</td><td>Bevacizumab, carboplatin, paclitaxel</td></tr><tr><td>NCT01004250</td><td>2</td><td>109</td><td>Bevacizumab, cisplatin, pemetrexed</td></tr><tr><td>NCT00312728</td><td>2</td><td>115</td><td>Bevacizumab, carboplatin, cisplatin, docetaxel, erlotinib, gemcitabine, paclitaxel, pemetrexed, vinorelbine</td></tr><tr><td>NCT00480831</td><td>2</td><td>128</td><td>Bevacizumab, carboplatin, paclitaxel, PRO95780 (apomab)</td></tr><tr><td>NCT00700180</td><td>2</td><td>303</td><td>Bevacizumab, carboplatin</td></tr><tr><td>NCT01364012</td><td>3</td><td>276</td><td>Bevacizumab, carboplatin, paclitaxel</td></tr><tr><td>NCT00762034</td><td>3</td><td>939</td><td>Bevacizumab, carboplatin, paclitaxel, pemetrexed</td></tr><tr><td>NCT00806923</td><td>3</td><td>1044</td><td>Bevacizumab, cisplatin, gemcitabine</td></tr><tr><td>NCT02596958</td><td>4</td><td>996</td><td>Bevacizumab, various platinum-based chemotherapeutics</td></tr><tr><td>NCT00388206</td><td>4</td><td>3998</td><td>Bevacizumab, various chemotherapeutics</td></tr><tr><td colspan="4"><strong>Additional Vivli studies not available in time for modeling</strong></td></tr><tr><td>doi:10.7150/jca.37966</td><td>N/A</td><td>55</td><td>Pembrolizumab, bevacizumab, nivolumab, various chemotherapeutics</td></tr><tr><td>NCT01846416</td><td>2</td><td>138</td><td>Atezolizumab</td></tr><tr><td>NCT01903993</td><td>2</td><td>287</td><td>Atezolizumab, docetaxel</td></tr><tr><td>NCT01984242</td><td>2</td><td>305</td><td>Atezolizumab, bevacizumab, sunitinib</td></tr><tr><td>NCT02031458</td><td>2</td><td>667</td><td>Atezolizumab</td></tr><tr><td>NCT02008227</td><td>3</td><td>1225</td><td>Atezolizumab, docetaxel</td></tr></tbody></table>

Note: Data was originally obtained from clinical study data research.com, but later during the project Vivli became the custodian of the datasets. During this transition, six more clinical trial datasets representing immune checkpoint inhibitors were added to the repository, but data were not able to be made anonymized and available in time for the modeling project. By following the links to the readers can learn further information about dataset demographics.
Due to contractual obligations between the data administrators (Roche, Eli Lilly) and data access provider (CSDR), the project files were moved to a new secure server managed by Vivli on September 22nd, 2020. During this transfer process, access to five additional datasets was requested involving newer standard-of-care therapeutics and NSCLC, but these data were not included in the final modeling dataset.

## 2.2 | Model Building

The key focal points of this project encompassed the precise modeling of individual tumor growth and response, the formulation of a mathematical framework that accounts for both acquired and intrinsic resistance, the characterization of individual drug effects and interactions, and the modeling of the temporal dynamics underlying tumor vascularization improvement through bevacizumab dosing.

Candidate models were built in multiple development phases. First, tumor PK profiles for each patient were generated using PK models and parameters from the literature. Model files were written using the Mlxtran language (Monolix Suite version 2020R1, Lixoft) [22]. Then, several sample datasets were created, with between 5% and 10% of the complete dataset for initial model building. This shortened calculation time and reduced computational complexity during model building. Next, several base candidate models were chosen for NSCLC pharmacodynamics (PD) and were finally implemented using the Mlxtran language.

Two well-established models of NSCLC were used as the base candidate structure of tumor growth and drug PD. The first model fit to the data was the Claret model of tumor growth and response to antiproliferatives [23]. The second model fit was built on a previous Gompertzian model of BEV-PEM/CIS [20, 24]. Initially, all antiproliferatives were assumed to operate in a simple log-kill manner [25].

Finally, using a combination of manual exploration and modified SAEM searches with simulated annealing on the $η_i$ terms, several reasonable parameter estimates with which to initialize the parameter search were identified. Lastly, the best-performing candidate models were finalized and compared for performance using quality-of-fit and robustness of parameter estimates, as previously described [26–29]. Any modifications to structure would be iterated on using a similar method and compared using a suite of model evaluation techniques.

## 2.3 | Clinical Trial Simulations

The objective of the clinical trial simulations was to assess the potential advantages of sequential administration of pemetrexed-cisplatin therapy with bevacizumab. For this purpose, a cohort of 5000 virtual individuals was simulated based on the final fit of the model using random draws from the individual parameter distributions. This cohort was subsequently replicated across a varying gap for sequential administration of bevacizumab and chemotherapy. The administration of bevacizumab followed the administration of pemetrexed and cisplatin with a time gap ranging from 0 to 3 days, in increments of 0.1 days (e.g., 0, 0.1,

0.2, ..., 2.8, 2.9, 3.0) resulting in the 31 sets of test conditions. These test conditions were derived from modifications of the Multicenter Phase II MAP Study of NSCLC [22]:

* Pemetrexed was administered every 22 days for 7 cycles at 945 mg and with an infusion time of 10 min (an average between the administration schedules in the datasets, although not aligned with the typical standard of care).

* Cisplatin was administered every 22 days for 3 cycles at 135 mg and with an infusion time of 30 min.

* Bevacizumab was administered every 22 days for 7 cycles at 795 mg and with an infusion time of 30 min.

This approach aimed to evaluate whether concurrent administration yielded superior outcomes compared to administering bevacizumab after a delayed interval following pemetrexed and cisplatin. All simulations were executed using Simulx (Monolix Suite version 2023R1, Lixoft). Subsequently, the most effective time gap was identified among all the gaps tested within the range of 0 to 3 days.

# 3 | Results

## 3.1 | Data Summary

Data were received in directories of fully anonymized data frames (e.g., Excel files, SAS files, *.csv, etc.) organized by general category of data. For initial data processing, data were explored for bulk trends, standardized, and data with suspected errors and outliers were removed. Studies were then collated in a single data frame designed for use with the Monolix suite (2020R1). Between the 11 studies, 3686 patients' data were determined to be potentially suitable for analysis. After exploring the datasets in greater detail, 2586 patients were determined to have all necessary data to create models of tumor growth and response, including unique patient IDs, a time recorded for each dosing and measurement event, individual tumor diameter measurements, and details of administration (Table 1). Due to the massive effort of this study and semi-mechanistic aims, demographic descriptions of the dataset are limited in this manuscript.

The longest diameter of each tumor was chosen as the dependent variable in the model. Statistically, all patient-tumor combinations were treated as unique subject-occasions. No distinction was made between within-subject variability and between-subject variability. While acknowledging the potential correlation of individual parameters across multiple tumors within an individual as well as variation across tumor location, it is important to note that the approach described above represents a compromise that strikes a favorable balance between modeling systems biology and model parameters structural identifiability—a chief example of the types of decisions made when constructing semi-mechanistic models. Among the subject-occasions, there were 6197 unique tumors belonging to 2586 patients.

After a short period of testing, several further restrictions were imposed on the dataset to facilitate numerical stability of the modeling. The first condition was that tumors were required to
have been measured three or more times to qualify for inclusion. This reduced the number of unique tumor IDs to 4701 and unique individuals to 2036. If for each sample $y_{ij}$, for individual $i$ at time $t$, $y_{ij}$ was greater than or equal to $y_{ij-1}$, the tumor was labeled as a monotonic non-responder and excluded from the dataset. Although monotonic non-responders may have been informative to the model, in practice, these IDs introduced computational errors into the Monolix routine. This reduced the number of individual tumors to 4473 and individual IDs to 1977. After removing these data, there were 4450 tumors and 1963 individuals left in the dataset. These restrictions imposed on the initial dataset representing 2586 patients reduced the number of samples from 29,885 to 26,515. This is an approximate 11% reduction in data. The reduction in the dataset size was primarily due to the limited frequency of tumor size assessments among the majority of patients in the clinical studies. Lastly, approximately 85% of the available data from each study (randomly allocated by subject-tumor pairs) was used for model building, and the remaining 15% of data was reserved for external validation of the final model.

## 3.2 | Model Building

Individual pharmacokinetic parameters were estimated in relatively few of the clinical trials referenced herein. Therefore, population pharmacokinetic models were collated from scientific literature (Table 2). Between-patient variability in drug pharmacokinetics was not included in the pharmacokinetic model as it led to practical unidentifiability of random effects parameters.

For pharmacodynamic modeling, two primary candidates were evaluated for the base description of tumor growth and response: (1) the Gompertzian model of tumor growth and (2) the Claret model of tumor growth [23, 25]. Both models describe tumor volumes, where the raw data reported the tumor's longest diameter. The tumor diameters were converted into volumes by making the assumption that the tumors were approximately spherical in

shape. Using the Bayesian information criteria as a parsimonious method of cross-evaluating models, the Gompertzian model outperformed the Claret model of tumor growth during several stages of structural model development.

In the Gompertzian model of tumor growth (Equation 1), the unperturbed tumor grows at a rate $\alpha$ and is exponentially limited in growth by parameter $\beta$. $v_c$ is a scaling factor relating individual tumor cell turnover to volume [24]. It was set to $10^6$ cells/mm$^3$, which is the classical assumption of the approximate number of cells per unit volume [14].

$$ \frac{dv}{dt} = \left( \alpha \cdot Q_{\alpha} - \beta \cdot \log\left(\frac{v}{v_c}\right) \right) \cdot v - Q_{\gamma} \cdot v $$ $$ \log(Q_{\alpha}) = -(1 + w_{\text{bev}_\alpha} \cdot \text{cc}_{\text{bev}}(t-\tau)) \cdot (\text{eff}_{\text{microt}} + \text{eff}_{\text{vegf}}) $$ $$ Q_{\gamma} = (1 + w_{\text{bev}_{\gamma}} \cdot \text{cc}_{\text{bev}}(t-\tau)) \cdot (\text{eff}_{\text{plat}} + \text{eff}_{\text{afolate}} + \text{eff}_{\text{dr5}} + \text{eff}_{\text{egfr}} + \text{eff}_{\text{dnasub}}) $$ (1)

$$ \frac{dv}{dt} = \left( \alpha - \beta \cdot \log\left(\frac{v}{v_0}\right) \right) \cdot v - Q_{\gamma} \cdot v - Q_{\rho} \cdot v + k k_{\text{inj}} \cdot v_{\text{inj}} \cdot p $$ $$ \frac{dv_{\text{inj}}}{dt} = Q_{\rho} \cdot v - k k_{\text{inj}} \cdot v_{\text{inj}} - Q_{\gamma} \cdot v_{\text{inj}} $$ $$ Q_{\rho} = (1 + w_{\text{bev}_{\rho}} \cdot \text{cc}_{\text{bev}}(t-\tau)) \cdot (\text{eff}_{\text{microt}} + \text{eff}_{\text{vegf}}) $$ $$ Q_{\gamma} = (1 + w_{\text{bev}_{\gamma}} \cdot \text{cc}_{\text{bev}}(t-\tau)) \cdot (\text{eff}_{\text{plat}} + \text{eff}_{\text{afolate}} + \text{eff}_{\text{dr5}} + \text{eff}_{\text{egfr}} + \text{eff}_{\text{dnasub}}) $$ (2)

$$ \frac{dz_1}{dt} = Q_{\gamma} \cdot v + Q_{\gamma} \cdot v_{\text{inj}} + k k_{\text{inj}} \cdot v_{\text{inj}} \cdot (1-p) - k k \cdot z_1 $$ $$ \frac{dz_2}{dt} = k k \cdot z_1 - k k \cdot z_2 $$ $$ \frac{dz_3}{dt} = k k \cdot z_2 - k k \cdot z_3 $$ $$ n = v + z_1 + z_2 + z_3 + v_l $$ $$ \text{Tumor Diameter} = 2 \cdot \left( \frac{3n}{4\pi} \right)^{\frac{1}{3}} $$ (3)

<table><caption>TABLE 2 | Pharmacokinetic parameters estimates.</caption><thead><tr><th>Therapeutic</th><th>V1</th><th>V2</th><th>V3</th><th>Q1</th><th>Q2</th><th>k12</th><th>k21</th><th>Cl</th><th>Source</th></tr></thead><tbody><tr><td>Units</td><td>L</td><td>L</td><td>L</td><td>L/day</td><td>L/day</td><td>day<sup>-1</sup></td><td>day<sup>-1</sup></td><td>L/day</td><td>—</td></tr><tr><td>Bevacizumab</td><td>2.80</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.223</td><td>0.215</td><td>0.216</td><td>[30]</td></tr><tr><td>Cisplatin</td><td>22.3</td><td>77.0</td><td>—</td><td>456.0</td><td>—</td><td>—</td><td>—</td><td>6.48</td><td>[31]</td></tr><tr><td>Pemetrexed</td><td>12.9</td><td>3.38</td><td>—</td><td>20.7</td><td>—</td><td>—</td><td>—</td><td>131.9</td><td>[32]</td></tr><tr><td>Apomab</td><td>3.97</td><td>3.84</td><td>—</td><td>0.793</td><td>—</td><td>—</td><td>—</td><td>0.328</td><td>[33]</td></tr><tr><td>Paclitaxel</td><td>229.0</td><td>856.0</td><td>30.3</td><td>3216.0</td><td>5112.0</td><td>—</td><td>—</td><td>10296.0</td><td>[34]</td></tr><tr><td>Carboplatin</td><td>11.9</td><td>8.23</td><td>—</td><td>2172.0</td><td>—</td><td>—</td><td>—</td><td>177.12</td><td>[35]</td></tr><tr><td>Gemcitabine</td><td>15.0</td><td>15.0</td><td>—</td><td>1008.0</td><td>—</td><td>—</td><td>—</td><td>3888.0</td><td>[36]</td></tr><tr><td>Docetaxel</td><td>7.9</td><td>—</td><td>—</td><td>—</td><td>—</td><td>27.12</td><td>3.6</td><td>723.12</td><td>[37]</td></tr><tr><td>Erlotinib<sup>a</sup></td><td>210.0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>102.96</td><td>[38, 39]</td></tr></tbody></table>

Note: V1 represents volume 1; V2 represents volume 2; V3 represents volume 3; Q1 represents intercompartmental clearance between volumes 1 and 2; Q2 represents intercompartmental clearance between volumes 2 and 3; k12 represents rate of transfer from volume 1 to volume 2; k21 represents rate of transfer from volume 2 to volume 1; Cl represents clearance from V1.

<footnote>aBioavailability estimated at 60%, ka estimated at 21.36 days<sup>-1</sup> [38, 39].</footnote>
In Equation (1), $Q_α$ and $Q_γ$ are the antiproliferative effects result in growth reduction and irreversible cell death, respectively. Chemotherapeutics, which acted on the microtubules—docetaxel and paclitaxel—along with the direct effect of bevacizumab, were included in $Q_α$. All other therapeutics were modeled as drugs resulting in irreversible cell death. Transient enhancement in efficacy via tumor vasculature normalization by bevacizumab was modeled as occurring at time $(t - τ)$ to account for the time delay between administration and efficacy enhancement, i.e., $τ$ [19].

Equation (1) heavily exaggerated the effect of bevacizumab in limiting cell growth rates. To account for this, a second compartment was implemented, which represented reversible cell injury. $v_{inj}$, from which cells could return to the primary volume (i.e., $v$) from (Equation 2). In this equation, the elimination rate from tumor cell injury is governed by the intercompartmental transfer rate $kk_{inj}$ as well as the proportion of repaired cells returned to the unperturbed cycle of proliferation, $p$ (value fixed between 0 and 1 using logit-link). Notably, in version 4b, the Gompertzian model growth term was scaled to the initial growth rate. Using log-rules, it can be shown that $log(V/V_c)—log(V_0/V_c)$ is equivalent to $log(V/V_0)$. This modified model structure provided a more accurate representation of the growth-limiting effects exerted by bevacizumab, paclitaxel, and docetaxel, striking a balance that avoided the exaggerated effect observed in Equation (1).

Lastly, resistance to anticancer treatment was modeled using a variant of an exposure-response model [40]. In this mathematical representation, the tumor cells were assumed to become increasingly more resistant to treatment (rate governed by parameter $λ_d$) as a function of drug exposure (area under the concentration time-curve, AUC) – Equation (6). As was done for the efficacy of the drugs, the effect of drugs with similar mechanisms of action was weighted relative to each other—i.e. $w_{d_{1-rel-2}}$ represents the weight of some drug (given the arbitrary index 1) relative to some other drug (drug 2) which shares a mechanism of action, e.g., cisplatin and carboplatin. This is meant to model the phenomenon whereby a tumor becomes resistant to a therapeutic mechanism of action rather than an individual drug. A complete diagram of the model is available for review in Figure 1.

$$ \begin{gathered} \text{eff}_d(t) = (w_d \cdot cc_d(t)) \\ \text{eff}_{\text{total}}(t) = \sum w_d \cdot cc_d(t) \end{gathered} \qquad (4) \\ d \in \text{cis, car, pem, apo, erl, gem, pac, doc, bev} $$

$$ \begin{gathered} \text{eff}_d(t) = (w_d \cdot cc_d(t)) \\ \text{eff}_{\text{total}}(t) = \gamma \cdot (\text{invlogit}(\sum w_d \cdot cc_d(t)) - 0.5) \end{gathered} \qquad (5) \\ d \in \text{cis, car, pem, apo, erl, gem, pac, doc, bev} $$

$$ \begin{gathered} \text{AUC}_d = \int cc_d(t) \\ \left[ \begin{array}{l} \text{singledrugeff}(t) \equiv \log(\text{eff}_m(t)) = \log(w_d \cdot cc_d(t)) - (\lambda_d \cdot \text{AUC}_d(t)) \\ \text{paireddrugeff}(t) \equiv \log(\text{eff}_m(t)) = \log(w_{d_{1-rel-2}}(cc_{d_1}(t) + w_{d_2} \cdot cc_{d_2}(t)) - \lambda_{d_{1-rel-2}}(\text{AUC}_{d_1}(t) + \lambda_{d_2} \cdot \text{AUC}_{d_2}(t)) \end{array} \right] \qquad (6) \\ d \in \text{cis, car, pem, apo, erl, gem, pac, doc, bev} \\ m \in \text{plat, afolate, dr5, egfr, dnasub, microt, vegf i. e. mode of action} \end{gathered} $$

Irreversible cell death was modeled as occurring over a series of transitions between several compartments with an intercompartmental transfer rate $kk$. The final tumor volume, a summation of the primary tumor volume and death compartments ($z_1$, $z_2$, and $z_3$), as well as the injured cell volume $v_{inj}$, was then transformed to a tumor diameter to match the (observed) independent variable in the dataset (Equation 3).

For the pharmacodynamic effect (eff) of the various therapeutics on tumor growth, a version of a previously published log-kill model was implemented, whereby each drug concentration was scaled by a single parameter (represented by $w_d$ for weighting for drug $d$) and the weighted sum of those concentrations determine the overall cytotoxic effect (Equation 4). $cc_d(t)$ represents the concentration of drug $d$ at any arbitrary time $t$. This formulation proved slightly unstable, so the sum of effects was eventually grouped under an inverse logit function, so that this sum would be bounded to a value between 0 and 1 (Equation 5). After testing several more versions of the model, it was found that grouping treatments together that shared a therapeutic mechanism-of-action (e.g., carboplatin and cisplatin) improved model stability and fit.

Individual variability was modeled using the standard lognormal distribution, and initial tumor volumes were fixed to the measurement of the tumor at time 0, relative to each individual. The only exception was made for parameter $p$, which was fixed between 0 and 1 using a logit-link function. Measurement error was modeled using the function *combined* 1 in Monolix, i.e., including a single additive term (a) added to a single proportional term (b).

## 3.3 | Model Diagnostics

Graphical evaluation of model diagnostics supports the validity and predictive performance of the final selected model. Specifically, the SAEM search was stable and reliable when estimating the final set of parameter estimates (see Figure S1). Individual fits were reasonably descriptive, with both a well-captured tumor growth inhibition in response to treatment and rebound after treatment cessation (see Figure S2). After evidence of correlations (via Pearson's test) between individual effects was observed, those correlations were investigated using the full posterior plot of individual effects. Although several potential correlations were discovered, the slopes of these correlations
<figure><img src="image_2.png" /><figcaption>**FIGURE 1** | Diagram of the systems pharmacology model. In the Gompertzian model of tumor growth, the unperturbed tumor grows at a rate $\alpha$ and is exponentially limited in growth by parameter $\beta$. When a cytotoxic is introduced into the system, the cytotoxic impairs the growth of the tumor by sending cells into a death succession governed by $kk$. The cytotoxic death relative to drug concentration is modeled as $Q_p$. Reversible cell injury relative to drug concentration is modeled as $Q_d$, $\tau$ is the time delay between bevacizumab administration and the perfusion enhancement effect in the tumor. When a cell is damaged by cytotoxics or reversible cell injury that progresses to irreversible cell injury, it begins a progression from unperturbed growth—compartment V(t)—to damage compartments Z1 through Z1. Eventually, the cell exits the tumor volume as it dies.</figcaption></figure>

were nearly zero, and they were likely detected in error (type 1) as an effect of working with such a large dataset (increased statistical power). An even spread of observations versus individual predictions suggests that this model has no major structural misspecifications and that the error model was well specified (see Figure S3). Formal tests for residual normality and centering on zero indicated that the model residual error matched theoretical expectations, other than predictions at time 0 and BLQ data. The non-normal distribution of residuals at time 0 and for BLQ data is due to the use of the initial measured tumor volume as the predicted tumor volume at time 0 (i.e., residual is equal to 0) and the conservative estimate of a lower limit of detection of 1 mm for measuring tumor diameter. The spread of individual parameters met the Kolmogorov-Smirnov test for normality. The precision of parameter estimates was extremely high (RSE < 8%) with low correlation between estimates. Shrinkage of individual parameter distributions toward the mean was low. Full parameter estimates, IIV, and RSEs are reported in Table 3. The visual predictive check (VPC) was informative as to the logistic model performance. Although the clinical trials were not matched in terms of sampling schedule, the VPC still indicates an overall high-quality fit (see Figure S4).

## 3.4 | External Validation

The model was externally validated using the reserved 15% of the data not used during the model-building process. Individual parameters were found for each subject, but population parameter estimates, IIV, and error parameters were all fixed to the value found in the final model fit. Model diagnostics were used to assess the ability of the model to make predictions on new patients that were not part of the original dataset. Using the VPC as an overall diagnostic tool (Figure 2), the final model fit the external data extremely well with little misspecification except in the

lower prediction interval band. This VPC was produced without re-estimating the population parameters, so it is a diagnostic for the fit of the model to an external dataset. Then, the reserved 15% of data was reduced to three subsets, one consisting of only the first datapoint for each patient, another consisting of only the first two datapoints, and a third consisting of only the first three datapoints. Individual predictions for each simulation scenario were compared between the three datasets to assess the quality of the model's individual predictions (Figure 3). The model could be used to predict the efficacy of individual treatment time courses, with an average difference between final prediction and observation of 59.7% after a single tumor measurement and 11.7% after three successive tumor measurements.

## 3.5 | Clinical Trial Simulations

Clinical trial simulations indicated that of the patients who did not go into complete remission in symptoms from standard concomitant therapy, a majority would have seen some improvement in tumor size reduction if they had received bevacizumab before pemetrexed-cisplatin. The ideal predicted gap, in steps of 0.1 days, for patients receiving bevacizumab-pemetrexed-cisplatin therapy was 0.4 days or 9.6 h (Figure 4). At this gap, approximately 93.5% of patients benefited from a gap in sequential administration. Of those patients, the mean improvement was 20.7%. The mean improvement in patients at the more practical gaps of 12 and 24 h was 15.7% and 14.3%, respectively.

# 4 | Discussion

The primary objective of this project was to develop a comprehensive semi-mechanistic model that encompasses the growth and response of non-small cell lung cancer (NSCLC) to various
<table><caption>TABLE 3 | Pharmacodynamic parameters estimates from the final model.</caption><thead><tr><th>Fixed effects</th><th>Value</th><th>SE</th><th>RSE (%)</th><th>Unit</th><th>IIV</th></tr></thead><tbody><tr><td>α</td><td>0.00304</td><td>8.27E-05</td><td>2.72</td><td>day<sup>-1</sup></td><td>1.22</td></tr><tr><td>β</td><td>0.0398</td><td>0.000546</td><td>1.37</td><td>day<sup>-1</sup></td><td>0.62</td></tr><tr><td>τ</td><td>0.254</td><td>0.00183</td><td>0.723</td><td>day</td><td>0.304</td></tr><tr><td>kk</td><td>28</td><td>2.05</td><td>7.3</td><td>day<sup>-1</sup></td><td>3.41</td></tr><tr><td>w<sub>cisplatin relative to carboplatin</sub></td><td>0.917</td><td>0.0728</td><td>7.94</td><td>—</td><td>3.59</td></tr><tr><td>w<sub>carboplatin</sub></td><td>0.257</td><td>0.00486</td><td>1.89</td><td>—</td><td>0.808</td></tr><tr><td>w<sub>pemetrexed</sub></td><td>23.8</td><td>0.759</td><td>3.19</td><td>—</td><td>1.41</td></tr><tr><td>w<sub>apomab</sub></td><td>0.75</td><td>0.00342</td><td>0.456</td><td>—</td><td>0.194</td></tr><tr><td>w<sub>erlotinib</sub></td><td>0.62</td><td>0.006</td><td>0.968</td><td>—</td><td>0.407</td></tr><tr><td>w<sub>gemcitabine</sub></td><td>1</td><td>0.0146</td><td>1.45</td><td>—</td><td>0.616</td></tr><tr><td>w<sub>paclitaxel relative to docetaxel</sub></td><td>0.758</td><td>0.00956</td><td>1.26</td><td>—</td><td>0.537</td></tr><tr><td>w<sub>docetaxel</sub></td><td>0.939</td><td>0.00233</td><td>0.248</td><td>—</td><td>0.105</td></tr><tr><td>w<sub>bevacizumab,ρ</sub></td><td>4.14</td><td>0.283</td><td>6.82</td><td>—</td><td>3.44</td></tr><tr><td>λ<sub>cisplatin relative to carboplatin</sub></td><td>0.00101</td><td>9.18E-06</td><td>0.907</td><td>—</td><td>0.386</td></tr><tr><td>λ<sub>carboplatin</sub></td><td>0.000995</td><td>9.58E-06</td><td>0.963</td><td>—</td><td>0.409</td></tr><tr><td>λ<sub>pemetrexed</sub></td><td>0.000796</td><td>7.25E-06</td><td>0.911</td><td>—</td><td>0.385</td></tr><tr><td>λ<sub>apomab</sub></td><td>0.000656</td><td>1.35E-06</td><td>0.206</td><td>—</td><td>0.0877</td></tr><tr><td>λ<sub>erlotinib</sub></td><td>0.00119</td><td>2.15E-05</td><td>1.81</td><td>—</td><td>0.777</td></tr><tr><td>λ<sub>gemcitabine</sub></td><td>0.00152</td><td>1.66E-05</td><td>1.09</td><td>—</td><td>0.47</td></tr><tr><td>λ<sub>paclitaxel relative to docetaxel</sub></td><td>0.00111</td><td>1.02E-05</td><td>0.918</td><td>—</td><td>0.386</td></tr><tr><td>λ<sub>docetaxel</sub></td><td>0.00134</td><td>1.04E-05</td><td>0.776</td><td>—</td><td>0.327</td></tr><tr><td>λ<sub>bevacizumab,ρ</sub></td><td>0.0164</td><td>0.000619</td><td>3.79</td><td>—</td><td>1.69</td></tr><tr><td>w<sub>bevacizumab,γ</sub></td><td>0.727</td><td>0.00638</td><td>0.878</td><td>—</td><td>0.373</td></tr><tr><td>w<sub>bevacizumab,ρ</sub></td><td>1.39</td><td>0.031</td><td>2.23</td><td>—</td><td>0.957</td></tr><tr><td>kk<sub>inj</sub></td><td>0.0746</td><td>0.00364</td><td>4.89</td><td>—</td><td>2.31</td></tr><tr><td>p</td><td>0.364</td><td>0.0072</td><td>1.98</td><td>day<sup>-1</sup></td><td>1.37</td></tr><tr><td>a</td><td>0.0457</td><td>0.00116</td><td>2.53</td><td>—</td><td>—</td></tr><tr><td>b</td><td>0.105</td><td>0.00105</td><td>1</td><td>—</td><td>—</td></tr></tbody></table>

Note: SE is the standard error of the estimate and RSE is the relative standard error of the estimate. IIV is the interindividual variability, or in practical terms the standard deviation of the random effects. $w_{drug}$ parameters are unitless and represent the weight of effect of the treatment modality on Gompertzian growth. $λ_d$ parameters are also unitless and represent the effect of AUC on the effectiveness of the relative treatment modalities. a and b are error model parameters and therefore have no IIV.

clinical antiproliferative treatments. The model produced in this study effectively incorporates the antiproliferative effects of 11 different anti-NSCLC therapeutics used in 11 clinical trials, accounting for both intrinsic and acquired resistance to anticancer therapy. To accomplish this, the most reliable published pharmacokinetic (PK) models for these therapeutics have been integrated into this model, opting not to estimate interindividual variability for PK. The final selected model comprehensively demonstrates transient enhancement of perfusion resulting from anti-VEGF therapy, particularly bevacizumab. The individual predictions derived from this model exhibit a relatively high degree of precision, effectively

capturing the well-established rebound in growth that follows the cessation of treatment in non-small cell lung cancer.

While modeling resistance to treatment, AUC was used as a surrogate of drug exposure to model acquired resistance. AUC serves as a meaningful metric for modeling drug exposure, as it is an easily measurable parameter in a clinical setting. It provides a comprehensive representation of the concentration-time profile of a drug, taking into account both the peak concentration and the duration of exposure. By quantifying the total drug exposure over a specified period, AUC enables comparisons across different dosing regimens and facilitates the evaluation
<figure><img src="image_3.png" /><figcaption>**FIGURE 2** | Visual predictive check using 15% of the full dataset (Reserved set: external validation). The empirical median tumor longest diameter (central blue line) and upper and lower empirical 90% quantiles (upper and lower blue lines) are predicted precisely by the model fit prediction intervals. The 90% prediction intervals for the median tumor longest diameter (blue band) and upper and lower 90% quantiles (red bands) capture the majority of the variation in the dataset.</figcaption></figure>

of drug efficacy and safety. Therefore, incorporating AUC into the modeling framework allowed for a robust characterization of drug exposure and its relationship to clinical outcomes. Intrinsic resistance has been folded into the distribution of $w_d$ (i.e., weighting) terms in this system of equations. A known weakness in this approach to modeling resistance is that, while capturing a distribution of resistance time courses with respect to exposure, this model has not accounted for the individual systems and pathophysiology that might better explain and characterize that resistance [41].

During this investigation, it was discovered that incorporating a second cytotoxic effect to account for reversible cellular injury proved essential in capturing both the direct effects of bevacizumab and the impact of medications known to induce reversible cell injury, such as paclitaxel and docetaxel. By employing a predominantly linear combination of differential equation terms, this model offers a seamless avenue for integrating additional drug candidates and further system-level details. Moreover, a set of parameter estimates has been meticulously derived and rigorously validated, and the model code has been made openly accessible (see data sharing statement) for future research endeavors focused on non-small cell lung cancer. This transparency and availability aim to facilitate scientific progress and enhance the understanding of NSCLC dynamics and therapeutic interventions.

Direct comparison between cancer model parameterization is not always possible, but the Gompertz model is relatively well-published in NSCLC, so growth and decay parameters can be compared. The growth parameter ($\alpha$) derived in this study is similar to those derived in other studies [42]. However, many

studies use the assumption that the limiting population of cells in NSCLC is 1E12, which would translate to a diameter of roughly 124 mm if the tumor were roughly spherical [43]. In this study, it was determined that the corresponding limiting diameter, given typical parameter values, would be 21 mm. This indicates that previous estimates of the typical limiting population of cells were too high and that a better estimate might be on the order of 5E9 cells. However, fixing initial tumor size to observed size may have caused an underestimation of true $\alpha$ and $\beta$ IIV. This was a necessary simplification to stabilize the estimation of $\alpha$ and $\beta$. These parameters are typically correlated in NLME Gompertzian models.

When evaluating the predictive performance of the model, strong evidence was observed in support of the model structure choice. Parameter estimates and individual predictions were made with high precision, which is likely due to the size of the dataset included. Model structure was based on biological mechanisms, and interpretation of parameters is relatively straightforward—e.g., $\lambda$ parameters define the rate of acquired resistance versus exposure. Mechanisms grounded in biology afford the model additional longevity. On parameter estimate interpretation, relatively simple naming heuristics were used to aid in interpretation. $w_d$ weight parameters drug action against the tumor, i.e., the larger the $w_d$ parameter, the larger the action the drug takes proportionally to both the tumor size and concentration of the drug in plasma. The parameter $p$ indicates the proportion of cells in the injured volume that return to unperturbed tumor growth. Unperturbed tumor growth is governed by parameters $\alpha$ (exponential rate of growth of tumor) and $\beta$ (exponential rate of decrease in growth rate due to nutrient supply limitations).
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** | Individual bayesian predictions after limited sampling. This test was used to assess overfitting in the model. The reserved 15% of data was reduced to three subsets, one consisting of only the first datapoint for each patient, another consisting of only the first two timepoints, and a third consisting of only the first three timepoints. Individual predictions for each were compared between the three datasets to assess the quality of the model's individual predictions. The model could be used to predict the efficacy of individual treatment time courses, with an average difference between final prediction and observation of 59.7% after a single measurement and 11.7% after three successive measurements.</figcaption></figure>

There are several weaknesses that need to be acknowledged in this approach and addressed in future studies. Large amounts of data on cutting-edge first-line therapies–i.e., pembrolizumab and other immune-checkpoint inhibitors–were not included at this time. There are also several covariates, such as weight, body surface area, gender, and concomitant disease, which have not been included or tested in the model. This modeling project was semi-mechanistic in its aims. This is partially due to the nature of measurements of tumor size made in a clinical trial. Those measurements are sparsely made and strongly influenced by the technicians performing the tumor measurements. There was a richness in sampling between clinical trials, but within any patient, sampling was infrequent. This resulted in imprecise modeling heuristics, such as modeling resistance as a function of AUC. However, the platform produced in this study is heavily informed by clinical practice and will make an impactful platform for further iteration informed by either more mechanistic mechanisms or AI modeling methods [44–47]. Lastly, there is some bias in the dataset design of the study. In the modeling dataset, there was moderate dropout due to patient death, patient withdrawal, occasionally overly sparse sampling schedules, and the infrequent occurrence of monotonic nonresponse. Although no data imputations (e.g., last observation carried forward method) were used, these realities likely marginally biased the parameter estimation toward overestimating the

efficacy of the therapeutics [48]. And because these sparse samples are being taken over many different experimental designs and a large set of IIV parameters, there is likely a mild IIV inflation despite acceptable individual fitting (see Figure S4).

One of the primary features that this model was designed to capture is the transient enhancement of drug delivery after bevacizumab administration. Theoretically, this transient enhancement drives the synergism between bevacizumab and other antiproliferatives. Potentially, some other biological phenomenon may drive the synergism between bevacizumab and other antiproliferatives, but this model is ultimately agnostic about the precise biological mechanism of synergism. A natural conclusion, and a finding supported by previously published clinical papers [19, 49–53], is that administering bevacizumab before other antiproliferatives should result in the greatest reduction in tumor size. A delay of 9.6 h between pemetrexed-cisplatin and bevacizumab resulted in the greatest benefit to the virtual patients. At this gap, approximately 93.5% of patients benefited from a gap in sequential administration. Of those patients, the mean improvement was 20.7%. The mean improvement in patients at the more practical gaps of 12 and 24 h was 15.7% and 14.3%, respectively. This is within a factor of 3 of the prediction of optimal gap made in Schneider et al. 2019 (1.2 days) using
<figure><img src="image_5.png" /><figcaption>**FIGURE 4** | Impact of gap in scheduling on overall efficacy. This stacked bar chart is meant to help determine the gap where simulated patients saw the best improvement from delayed bevacizumab treatment. The shaded regions represent the percentage of simulated patients who saw improvement over concomitant dosing, quantified by reduction in final tumor size, for the specified gap in days. The highest intensity of improvement is at a gap of 9.6 h, or 0.4 days. On the left-hand side, the axis represents the percentage of patients. After approximately 5 days, there is no longer any predicted improvement gained from delaying bevacizumab administration.</figcaption></figure>

experimental data from xenograft mice [20]. In the future, this model could be adapted to meet standard clinical endpoints of late-stage NSCLC clinical trials, such as survival (OS) and progression-free survival (PFS), with RECIST criteria [54].

## 5 | Conclusion and Future Directions

This research provides a simulation platform for future *in silico* studies of optimal scheduling of various therapeutics. The individual Bayesian predictions made in the reserved 15% of data showed that the model had little overfitting and can be extended to external datasets. Mean errors in individual predictions with 1, 2, or 3 sampled points were extremely low, even for far-off successive timepoints. For future clinical therapy, this study reaffirms, in human patients, the advantage of staggering bevacizumab, pemetrexed, and cisplatin therapy. Lastly, this platform is well supported within the context of this study, but with careful refinement and further validation, this model is flexible enough to be adapted to accelerate preclinical predictions in related antiproliferatives. Another natural extension would be to explicitly model the variability due to tumor location as well as within-subject variability. Further research using newer modalities such as immune checkpoint inhibitors and newer methods such as AI-powered modeling will benefit greatly from this significant modeling effort.

### Author Contributions

All authors wrote the manuscript and designed the research. B.K.S. performed the research and analyzed the data. PI/Committee, review, oversight: Jonathan P. Mochel and Sebastien Benzekry. All other research was a joint effort with B.K.S. as the primary researcher.

### Acknowledgments

This manuscript is based on research using data from data contributors Roche and Eli Lilly that has been made available through Vivli Inc. Vivli has not contributed to or approved, and Vivli, Eli Lilly, and Roche are not in any way responsible for the contents of this publication.

### Conflicts of Interest

The authors declare no conflicts of interest.

### Data Availability Statement

The data that support the findings of this study are available from Vivli Center for Global Clinical Research Data, but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Some simulated data are, however, available from the authors upon reasonable request to Benjamin K. Schneider and with permission of the Vivli Center for Global Clinical Research Data.

### References

1. J. Ferlay, M. Ervik, F. Lam, et al., "Global Cancer Observatory: Cancer Today," Lyon, France: International Agency for Research on Cancer, accessed November 20, 2020. Cancer Today, http://gco.iarc.fr/today/home.

2. V. T. DeVita, T. S. Lawrence, and S. A. Rosenberg, eds., *DeVita, Hellman, and Rosenberg's Cancer: Principles & Practice of Oncology*, 11th ed. (Wolters Kluwer, 2019), 2390.

3. J. E. Niederhuber, *Abeloff's Clinical Oncology*, 6th ed. (Elsevier, 2019).

4. National Cancer Institute (NCI): Surveillance, Epidemiology, and End Results Program (SEER), "SEER*Explorer: An Interactive Website for SEER Cancer Statistics," accessed February 5, 2021, https://seer.cancer.gov/explorer/index.html.
5. American Cancer Society, "Cancer Facts & Statistics," accessed February 14, 2021, http://cancerstatisticscenter.cancer.org/.

6. D. S. Ettinger, D. E. Wood, D. L. Aisner, et al., "National Comprehensive Cancer Network Guidelines in Oncology: Non-Small Cell Lung Cancer Guidelines With NCCN Evidence Blocks (Version 2.2021—January 5, 2021)," National Comprehensive Cancer Network (2021), https://www.nccn.org/professionals/physician_gls/pdf/nscl_blocks.pdf.

7. N. H. Hanna, B. J. Schneider, S. Temin, et al., "Therapy for Stage IV Non-Small-Cell Lung Cancer Without Driver Alterations: ASCO and OH (CCO) Joint Guideline Update," *Journal of Clinical Oncology* 38, no. 14 (2020): 1608-1632, https://doi.org/10.1200/JCO.19.03022.

8. G. A. Masters, S. Temin, C. G. Azzoli, et al., "Systemic Therapy for Stage IV Non-Small-Cell Lung Cancer: American Society of Clinical Oncology Clinical Practice Guideline Update," *Journal of Clinical Oncology* 33, no. 30 (2015): 3488-3515.

9. M. A. Socinski, R. M. Jotte, F. Cappuzzo, et al., "Atezolizumab for First-Line Treatment of Metastatic Nonsquamous NSCLC," *New England Journal of Medicine* 378, no. 24 (2018): 2288-2301.

10. L. Gandhi, D. Rodríguez-Abreu, S. Gadgeel, et al., "Pembrolizumab Plus Chemotherapy in Metastatic Non-Small-Cell Lung Cancer," *New England Journal of Medicine* 378, no. 22 (2018): 2078-2092.

11. S. Gadgeel, D. Rodríguez-Abreu, G. Speranza, et al., "Updated Analysis From KEYNOTE-189: Pembrolizumab or Placebo Plus Pemetrexed and Platinum for Previously Untreated Metastatic Nonsquamous Non-Small-Cell Lung Cancer," *Journal of Clinical Oncology* 38, no. 14 (2020): 1505-1517.

12. T. A. d'Amato, R. J. Landreneau, R. J. McKenna, R. S. Santos, and R. J. Parker, "Prevalence of In Vitro Extreme Chemotherapy Resistance in Resected Nonsmall-Cell Lung Cancer," *Annals of Thoracic Surgery* 81, no. 2 (2006): 440-447.

13. W. Pao, T. Y. Wang, G. J. Riely, et al., "KRAS Mutations and Primary Resistance of Lung Adenocarcinomas to Gefitinib or Erlotinib," *PLoS Medicine* 2, no. 1 (2005): e17.

14. Z. H. Tang, X. M. Jiang, X. Guo, C. M. V. Fong, X. Chen, and J. J. Lu, "Characterization of Osimertinib (AZD9291)-Resistant Non-Small Cell Lung Cancer NCI-H1975/OSIR Cell Line," *Oncotarget* 7, no. 49 (2016): 81598-81610.

15. E. B. Garon, M. D. Hellmann, N. A. Rizvi, et al., "Five-Year Overall Survival for Patients With Advanced Non-Small-Cell Lung Cancer Treated With Pembrolizumab: Results From the Phase I KEYNOTE-001 Study," *Journal of Clinical Oncology* 37, no. 28 (2019): 2518-2527.

16. T. J. Lynch, D. R. Spigel, J. Brahmer, et al., "Safety and Effectiveness of Bevacizumab-Containing Treatment for Non-Small-Cell Lung Cancer: Final Results of the ARIES Observational Cohort Study," *Journal of Thoracic Oncology* 9, no. 9 (2014): 1332-1339.

17. A. M. Stein and M. Looby, "Benchmarking QSP Models Against Simple Models: A Path to Improved Comprehension and Predictive Performance," *CPT: Pharmacometrics & Systems Pharmacology* 7, no. 8 (2018): 487-489.

18. A. Joshi, S. Ramanujan, and J. Y. Jin, "The Convergence of Pharmacometrics and Quantitative Systems Pharmacology in Pharmaceutical Research and Development," *European Journal of Pharmaceutical Sciences* 182 (2023): 106380.

19. D. C. Imbs, R. E. Cheikh, A. Boyer, et al., "Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma," *CPT: Pharmacometrics & Systems Pharmacology* 7, no. 1 (2018): 42-50.

20. B. Schneider, B. Bieth, A. Boyer, et al., "Optimal Scheduling of Bevacizumab and Pemetrexed/Cisplatin Dosing in Non-Small Cell Lung Cancer," *Abstracts of the Annual Meeting of the Population Approach Group in Europe* 28 (2019).

21. Vivli, "Vivli—Center for Global Clinical Research Data," accessed February 7, 2021, https://vivli.org/.

22. "Mlxtran|A Human Readable Language for Model Description," accessed April 23, 2023, https://mlxtran.lixoft.com/.

23. L. Claret, P. Girard, P. M. Hoff, et al., "Model-Based Prediction of Phase III Overall Survival in Colorectal Cancer on the Basis of Phase II Tumor Dynamics," *Journal of Clinical Oncology* 27, no. 25 (2009): 4103-4108.

24. C. Vaghi, A. Rodallec, R. Fanciullino, et al., "Population Modeling of Tumor Growth Curves and the Reduced Gompertz Model Improve Prediction of the Age of Experimental Tumors," *PLoS Computational Biology* 16, no. 2 (2020): e1007178.

25. S. Benzekry, E. Pasquier, D. Barbolosi, et al., "Metronomic Reloaded: Theoretical Models Bringing Chemotherapy Into the Era of Precision Medicine," *Seminars in Cancer Biology* 35 (2015): 53-61.

26. L. Pelligand, A. Soubret, J. N. King, J. Elliott, and J. P. Mochel, "Modeling of Large Pharmacokinetic Data Using Nonlinear Mixed-Effects: A Paradigm Shift in Veterinary Pharmacology. A Case Study With Rovenabexib in Cats," *CPT: Pharmacometrics & Systems Pharmacology* 5, no. 11 (2016): 625-635.

27. J. Wang, B. K. Schneider, J. Xue, et al., "Pharmacokinetic Modeling of Ceftiofur Sodium Using Nonlinear Mixed-Effects in Healthy Beagle Dogs," *Frontiers in Veterinary Science* 6 (2019): 363, https://doi.org/10.3389/fvets.2019.00363/full.

28. J. Wang, B. K. Schneider, H. Xiao, et al., "Non-Linear Mixed-Effects Pharmacokinetic Modeling of the Novel COX-2 Selective Inhibitor V泛acoxib in Cats," *Frontiers in Veterinary Science* 7 (2020): 554033.

29. J. Wang, B. K. Schneider, P. Sun, et al., "Nonlinear Mixed-Effects Pharmacokinetic Modeling of the Novel COX-2 Selective Inhibitor V泛acoxib in Dogs," *Journal of Veterinary Pharmacology and Therapeutics* 42, no. 5 (2019): 530-540.

30. J. F. Lu, R. Bruno, S. Eppler, W. Novotny, B. Lum, and J. Gaudreault, "Clinical Pharmacokinetics of Bevacizumab in Patients With Solid Tumors," *Cancer Chemotherapy and Pharmacology* 62, no. 5 (2008): 779-786.

31. S. Urien, E. Brain, R. Bugat, et al., "Pharmacokinetics of Platinum After Oral or Intravenous Cisplatin: A Phase 1 Study in 32 Adult Patients," *Cancer Chemotherapy and Pharmacology* 55, no. 1 (2005): 55-60.

32. J. E. Latz, A. Chaudhary, A. Ghosh, and R. D. Johnson, "Population Pharmacokinetic Analysis of Ten Phase II Clinical Trials of Pemetrexed in Cancer Patients," *Cancer Chemotherapy and Pharmacology* 57, no. 4 (2006): 401-411.

33. D. R. Camidge, R. S. Herbst, M. S. Gordon, et al., "A Phase I Safety and Pharmacokinetic Study of the Death Receptor 5 Agonistic Antibody PRO95780 in Patients With Advanced Malignancies," *Clinical Cancer Research* 16, no. 4 (2010): 1256-1263.

34. A. Henningsson, M. O. Karlsson, L. Viganò, L. Gianni, J. Verweij, and A. Sparreboom, "Mechanism-Based Pharmacokinetic Model for Paclitaxel," *Journal of Clinical Oncology* 19, no. 20 (2001): 4065-4073.

35. M. Joerger, A. D. R. Huitema, D. J. Richel, et al., "Population Pharmacokinetics and Pharmacodynamics of Paclitaxel and Carboplatin in Ovarian Cancer Patients: A Study by the European Organization for Research and Treatment of Cancer-Pharmacology and Molecular Mechanisms Group and New Drug Development Group," *Clinical Cancer Research* 13, no. 21 (2007): 6410-6418.

36. X. Jiang, P. Galettis, M. Links, P. L. Mitchell, and A. J. McLachlan, "Population Pharmacokinetics of Gemcitabine and Its Metabolite in Patients With Cancer: Effect of Oxaliplatin and Infusion Rate," *British Journal of Clinical Pharmacology* 65, no. 3 (2008): 326-333.

37. K. A. Slaviero, S. J. Clarke, A. J. McLachlan, E. Y. L. Blair, and L. P. Rivory, "Population Pharmacokinetics of Weekly Docetaxel in Patients With Advanced Cancer," *British Journal of Clinical Pharmacology* 57, no. 1 (2004): 44-53.

38. J. F. Lu, S. M. Eppler, J. Wolf, et al., "Clinical Pharmacokinetics of Erlotinib in Patients With Solid Tumors and Exposure-Safety
Relationship in Patients With Non-Small Cell Lung Cancer,” *Clinical Pharmacology and Therapeutics* 80, no. 2 (2006): 136–145.

39. P. Frohna, J. Lu, S. Eppler, et al., “Evaluation of the Absolute Oral Bioavailability and Bioequivalence of Erlotinib, an Inhibitor of the Epidermal Growth Factor Receptor Tyrosine Kinase, in a Randomized, Crossover Study in Healthy Subjects,” *Journal of Clinical Pharmacology* 46, no. 3 (2006): 282–290.

40. R. V. Overgaard, S. H. Ingwersen, and C. W. Tornøe, “Establishing Good Practices for Exposure-Response Analysis of Clinical Endpoints in Drug Development,” *CPT: Pharmacometrics & Systems Pharmacology* 4, no. 10 (2015): 565–575.

41. M. Nagase, S. Aksenov, H. Yan, J. Dunyak, and N. Al-Huniti, “Modeling Tumor Growth and Treatment Resistance Dynamics Characterizes Different Response to Gefitinib or Chemotherapy in Non-Small Cell Lung Cancer,” *CPT: Pharmacometrics & Systems Pharmacology* 9, no. 3 (2020): 143–152.

42. P. Schlicke, C. Kuttler, and C. Schumann, “How Mathematical Modeling Could Contribute to the Quantification of Metastatic Tumor Burden Under Therapy: Insights in Immunotherapeutic Treatment of Non-Small Cell Lung Cancer,” *Theoretical Biology & Medical Modelling* 18, no. 1 (2021): 11.

43. J. S. Spratt and J. A. Spratt, “The Prognostic Value of Measuring the Gross Linear Radial Growth of Pulmonary Metastases and Primary Pulmonary Cancers,” *Journal of Thoracic and Cardiovascular Surgery* 71, no. 2 (1976): 274–278.

44. B. Ocaña-Tienda, J. Pérez-Beteta, J. D. Villanueva-García, et al., “A Comprehensive Dataset of Annotated Brain Metastasis MR Images With Clinical and Radiomic Data,” *Scientific Data* 10, no. 1 (2023): 1–6.

45. S. Haupt, A. Zeilmann, A. Ahadova, et al., “Mathematical Modeling of Multiple Pathways in Colorectal Carcinogenesis Using Dynamical Systems With Kronecker Structure,” *PLoS Computational Biology* 17, no. 5 (2021): e1008970.

46. A. L. Jenner, T. Cassidy, K. Belaid, M. C. Bourgeois-Daigneault, and M. Craig, “In Silico Trials Predict That Combination Strategies for Enhancing Vesicular Stomatitis Oncolytic Virus Are Determined by Tumor Aggressivity,” *Journal for Immunotherapy of Cancer* 9, no. 2 (2021): e001387.

47. pyDarwin Certara USA Inc. (2023), https://github.com/certara/pyDarwin.

48. C. Hu and M. E. Sale, “A Joint Model for Nonlinear Longitudinal Data With Informative Dropout,” *Journal of Pharmacokinetics and Pharmacodynamics* 30, no. 1 (2003): 83–103.

49. P. V. Dickson, J. B. Hamner, T. L. Sims, et al., “Bevacizumab-Induced Transient Remodeling of the Vasculature in Neuroblastoma Xenografts Results in Improved Delivery and Efficacy of Systemically Administered Chemotherapy,” *Clinical Cancer Research* 13, no. 13 (2007): 3942–3950.

50. J. Ciccolini, S. Benzekry, B. Lacarelle, D. Barbolosi, and F. Barlési, “Improving Efficacy of the Combination Between Antiangiogenic and Chemotherapy: Time for Mathematical Modeling Support,” *Proceedings of the National Academy of Sciences of the United States of America* 112, no. 27 (2015): E3453.

51. S. Mollard, J. Ciccolini, D. C. Imbs, R. E. Cheikh, D. Barbolosi, and S. Benzekry, “Model Driven Optimization of Antiangiogenics + Cytotoxics Combination: Application to Breast Cancer Mice Treated With Bevacizumab + Paclitaxel Doublet Leads to Reduced Tumor Growth and Fewer Metastasis,” *Oncotarget* 8, no. 14 (2017): 23087–23098.

52. A. A. M. Van der Veldt, M. Lubberink, I. Bahce, et al., “Rapid Decrease in Delivery of Chemotherapy to Tumors After Anti-VEGF Therapy: Implications for Scheduling of Anti-Angiogenic Drugs,” *Cancer Cell* 21, no. 1 (2012): 82–91.

53. S. M. Tolaney, Y. Boucher, D. G. Duda, et al., “Role of Vascular Density and Normalization in Response to Neoadjuvant Bevacizumab and Chemotherapy in Breast Cancer Patients,” *PNAS* 112, no. 46 (2015): 14325–14330.

54. Y. Wang, C. Sung, C. Dartois, et al., “Elucidation of Relationship Between Tumor Size and Survival in Non-Small-Cell Lung Cancer Patients Can Aid Early Decision Making in Clinical Drug Development,” *Clinical Pharmacology and Therapeutics* 86, no. 2 (2009): 167–174.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section.