Mathematical modeling identifies optimum lapatinib dosing schedules for the treatment of glioblastoma patients

Human primary glioblastomas (GBM) often harbor mutations within the epidermal growth factor receptor (EGFR). Treatment of EGFR-mutant GBM cell lines with the EGFR/HER2 tyrosine kinase inhibitor lapatinib can effectively induce cell death in these models. However, EGFR inhibitors have shown little efficacy in the clinic, partly because of inappropriate dosing. Here, we developed a computational approach to model the in vitro cellular dynamics of the EGFR-mutant cell line SF268 in response to different lapatinib concentrations and dosing schedules. We then used this approach to identify an effective treatment strategy within the clinical toxicity limits of lapatinib, and developed a partial differential equation modeling approach to study the in vivo GBM treatment response by taking into account the heterogeneous and diffusive nature of the disease. Despite the inability of lapatinib to induce tumor regressions with a continuous daily schedule, our modeling approach consistently predicts that continuous dosing remains the best clinically feasible strategy for slowing down tumor growth and lowering overall tumor burden, compared to pulsatile schedules currently known to be tolerated, even when considering drug resistance, reduced lapatinib tumor concentrations due to the blood brain barrier, and the phenotypic switch from proliferative to migratory cell phenotypes that occurs in hypoxic microenvironments. Our mathematical modeling and statistical analysis platform provides a rational method for comparing treatment schedules in search for optimal dosing strategies for glioblastoma and other cancer types.

Introduction Glioblastoma is the most common and aggressive form of brain tumors in adults, characterized by short survival and poor treatment response [1]. Currently, the standard of care for glioblastoma patients includes surgery followed by radiotherapy and adjuvant chemotherapy with temozolomide [2]. However, the addition of chemotherapy only modestly prolongs survival (median 14.6 months) compared to radiation alone (median 12.1 months). Thus, there remains a pressing unmet medical need for more effective therapeutic agents. Unfortunately, since the introduction of temozolomide, no other compound has been able to significantly prolong patient survival in clinical trials. For orally administered drugs, most trials have only explored daily continuous dosing schedules (Table 1). However, there is increasing evidence that for some targeted agents, intermittent schedules can deliver equal or potentially even superior therapeutic benefit with less toxicity [3,4].
In the last decade, several molecularly targeted agents that inhibit recurrently mutated proteins have been investigated as a therapeutic strategy in glioblastoma. These have included several inhibitors of the epidermal growth factor receptor (EFGR), which is mutationally activated in approximately 50% of adult GBMs. Lapatinib is a small molecule tyrosine kinase inhibitor of human epidermal receptor 2 (HER2) and EGFR, which currently has regulatory approval for the treatment of HER2-positive advanced or metastatic breast cancer [5]. Additional indications for lapatinib [6][7][8][9][10][11][12][13][14] including glioblastoma, are currently being explored [15,16]. In GBM, clinical trials of lapatinib have failed to show efficacy using continuous dosing [15,16]. Interestingly, a study that evaluated EGFR inhibition as a means to prime tumor vasculature for efficient delivery of chemotherapy showed that in glioblastoma patients, a 2-day pulse of high dose (5,250 mg/day) lapatinib given through twice daily dosing was well  [17]. However, the question remains whether an altered dosing strategy might increase the efficacy of this agent. Due to ethical concerns of testing several dosing strategies in the absence of preclinical data suggesting their benefit, as well as to speed up discovery, mathematical modeling of treatment response can be used to identify predicted best administration schedules. Here, we explored continuous and pulsatile dosing strategies in an EGFR-mutant GBM cell line, and used mathematical and statistical modeling to determine optimal lapatinib dosing schedules for inhibiting tumor growth. The use of mathematical models for treatment optimization is part of a growing effort to improve clinical trial design for cancer patients. Many models have been developed to investigate the relationship between toxicity and dose in an attempt to better characterize the toxicity profile of the therapeutic products under investigation. For instance, Huitema et al. reviewed several differential equation-based models of chemotherapy-induced myelosuppression, cardiovascular events and other ordinal adverse events [18]. Amantea et al. developed a model incorporating exposure, biomarkers, efficacy endpoints and adverse events in models describing the treatment response of patients with gastrointestinal stromal tumors [19]. Their analysis uncovered a correlation between adverse events and efficacy: drug-induced increases in diastolic blood pressure were positively associated with overall survival. Another study by Huitema et al. examined the effects of the anti-angiogenic drug E7080 and its drug-induced adverse events in an attempt to evaluate dosing regimens with regard to their reduction of adverse events and improvement of efficacy [20]. They found that proteinuria could be best described by a discrete-time Markov transition model. Similarly, Fuhr et al. constructed a continuous-time Markov model to investigate erlotinib-induced adverse events in non-small cell lung cancer patients; their simulation results provided support for the use of high-dose pulses as an alternative dosing strategy for addressing acquired resistance often found when using a low continuous dose [21]. In these examples, a delicate balance between toxicity and efficacy determines the most efficacious dose.
In addition to understanding toxicity profiles, mathematical models have also been used to model pre-clinical data from animal or in vitro systems in order to bypass the difficulty associated with directly observing tumor progression in humans. Translational models incorporating in vitro or animal data have also be developed to investigate the relationship between efficacy endpoints and exposure. Mould et al. provided a comprehensive review of how the use of mathematical models can aid early development of anti-cancer therapy, in particular for describing tumor volume as a function of drug exposure [22]. Most papers these authors discussed describe changing tumor volumes using differential equation models with two terms: net tumor growth in the absence of therapeutic candidates and a drug-induced shrinkage effect [23][24][25][26][27]. These models suggest that the best strategy for inhibiting tumor growth is to maximize the effect of drug exposure, which often corresponds to higher drug concentrations. However, given toxicity constraints, high drug concentrations usually cannot be achieved for long time periods. An alternative strategy is to administer drugs in pulsatile doses in order to circumvent the toxicity limit. A small trial using weekly pulsatile high-dose erlotinib showed promising results for controlling central nervous system metastases from epidermal growth factor receptor-mutant non-small cell lung cancer [28]. Crooke et al. also performed computational analyses comparing continuous and pulsatile dosing and found that pulsed therapy is more effective than continuous therapy [29]. A clinical trial [30] based on mathematical modeling and preclinical [31] experiments demonstrated that a combined low-dose continuous and high-dose pulsed erlotinib schedule was successful at preventing progression in patients with CNS metastases but did not show significantly delayed emergence of resistance due to the T790M EGFR mutation.
Here we present a computational modeling approach incorporating both toxicity data from phase I clinical trials and efficacy data from pre-clinical research to examine the effects of various lapatinib dosing schedules on tumor growth (Fig 1). The goal of our work is to determine the best dosing strategy given the observed toxicity limit. Specifically, we are interested in answering the question whether alternative dosing strategies can be applied to circumvent the toxicity limit while exceeding the level of efficacy observed in continuous dosing of lapatinib while taking into account known characteristics of GBM growth and treatment response such as diffusivity, intratumor heterogeneity, and the blood brain barrier.

Cell culture
SF268 cells were obtained from the National Cancer Institute (NIH) and routinely grown in DMEM supplemented with 10% fetal bovine serum. To assess the effects of drug treatment, 150,000 cells were seeded on 6-cm dishes and allowed to attach overnight. Cells were then switched to growth media supplemented with 1% fetal bovine serum and the indicated concentrations of lapatinib or DMSO (vehicle). Each treatment condition was done in triplicate. For treatment discontinuation, cells were washed 3 times with media containing vehicle, and then allowed to continue growth in media containing vehicle. Vehicle-treated cells were also washed 3 times to control for the effects of washing. Cell viability and cell death were evaluated by the trypan blue exclusion assay using a ViCell cell viability analyzer (Beckman Coulter).

Mathematical modeling
In vitro cell dynamics. To describe in vitro cell dynamics, we designed a differential equation model of cell growth given by where Y 1 and Y 2 denote the numbers of viable and dead cells, respectively; the parameters b and d denote the birth and death rates of viable cells; K denotes the carrying capacity of the in vitro assay system; and c denotes the clearance rate of dead cells. which minimizes the squared distance between the observed and predicted cell numbers for each experimental setting in order to estimate model parameters that best recapitulate the experimental results. Characterizing the relationship between birth and death rates and lapatinib dose. We used an exponential model to describe the relationship between birth and death rates and lapatinib dose given by where b 0 and b 1 are the intercept and slope parameters, respectively, describing the birth rate, and d 0 and d 1 are the intercept and slope parameters, respectively, describing the death rate. Patient-derived toxicity constraints. To characterize the relationship between the maximally tolerated lapatinib dose per treatment and number of treatment days, we fit a linear model given by and an exponential model given by where y dose denotes the lapatinib intake per treatment day and x days is the number of treatment days per treatment cycle. Lapatinib pharmacokinetics. We fit a logarithmic function to model lapatinib absorption into the blood stream, according to where x oral is the dose of oral lapatinib in miligrams and y serum is the concentration of lapatinib in the serum in ng/mL units.
To convert serum to tumor lapatinib concentration, we used the conversion where y serum is the serum concentration in ng/mL, y tumor is the concentration of lapatinib in the tumor, and 0.61 was the mean tumor/serum ratio reported in Vivanco et al. [32].
Modeling in vivo GBM growth using a diffusion equation. We used the following diffusion model to simulate in vivo GBM growth where c(x, t) is the concentration of cells in 1mm 3 volume of tissue at position x and time t, D is the diffusion coefficient defining the net rate of migration of tumor cells, ρ is the net proliferation rate of cells per day, κ is the carrying capacity of the tissue, and r 2 is the Laplacian.
To model GBM growth during treatment, we scale ρ by the birth and death rates estimated from the in vitro lapatinib treatment data such that where ρ 0 = 0.0387 is the proliferation rate in the absence of treatment [33], L(t) is the concentration of lapatinib in the tumor at time t, b(0) and d(0) are the birth and death rates in the absence of treatment, and b(L(t)) and d(L(t)) are the birth and death rates during treatment with L(t) concentration of lapatinib. The diffusion model during treatment is then given by Lastly, we assumed an isolated system corresponding to a no-flux boundary condition of the form @cð0; tÞ @x ¼ @cðR; where T is the simulation time and R is the radius of the computational domain.
Modeling GBM heterogeneity using diffusion equations. We used the following system of diffusion equations to model intratumor heterogeneity due to the presence of drug-sensitive and -resistant cells: where c S and c R represent the concentration of sensitive and resistant cells at position x and time t. The proliferation rate in Eq (12b) given by ρ 0 , the proliferation rate in the absence of treatment, because the resistant cells are assumed to be unaffected by treatment.
Modeling reduced lapatinib delivery to the brain. To model reduced lapatinib delivery to the brain due to the blood brain barrier or other reasons, we used a modified diffusion model given by where Q denotes the percent of the serum concentration of lapatinib that reaches the tumor. Modeling the Go-or-Grow mechanism. It is well known that glioblastoma is highly hypoxic, and that one mechanism used by tumor cells to survive in hypoxic environments is a phenotypic switch from proliferative to migratory phenotypes [34], also known as the "Go-or-Grow" mechanism. This phenotypic switch depends on the oxygen concentration in the microenvironment, which we modeled using the following system of three PDEs describing migratory cell, proliferative cell, and oxygen concentration, repectively [35]: where c M (x, t) and c P (x, t) represent the concentration of migratory and proliferative tumor cells, respectively, at position x and time t, and σ(x, t) is the oxygen concentration at position x and time t. The latter is given by the diffusion equation where D σ is the oxygen diffusion coefficient, h 1 is the transvasculature permeability, h 2 is the rate at which glioma cells consume oxygen, and σ 0 is the oxygen concentration in normal brain cells. Lastly, f MP (σ) and f PM (σ) represent the rates at which tumor cells switch from migratory to proliferative and proliferative to migratory phenotypes, respectively. Assuming the rates of switching have a linear dependence on oxygen concentration, the two functions are then given by where λ 1 and λ 2 are positive constants.
Testing the relationship between cell motility and treatment response. To test whether the relationship between the cell motility predicted by the diffusion equation and the tumor volume depends on the treatment schedule, we used a linear regression model with interaction terms: Here D is the diffusion parameter, S is a categorical variable representing the treatment schedule used, and V S is the volume at the end of treatment with treatment schedule S. Then, (ν D + ν int, S S) is the slope representing the effect of motility on the tumor volume at the end of treatment, and ν int, S S is the effect of the interaction between cell motility and the treatment schedule. We used schedule 1 as the baseline schedule, meaning that ν S = ν int, S = 0. Thus we can test H 0 : ν int, S = 0 to assess whether tumor growth under schedule S is more or less affected by increased cell motility than schedule 1, for S = schedule 2, schedule 3, schedule 4, or schedule 5.
Simulating a clinical trial of treatment schedules. A cohort of N = 50 patients was used to simulate each clinical trial to test for differences between treatment schedules. Cell dynamics and pharmacokinetics parameters were randomly sampled for each patient. The same 50 patients were used to evaluate survival outcomes for different treatment schedules. Cell birth and death rate parameters, b 0 , b 1 , d 0 , and d 1 , were sampled from Normalðm;ŝÞ distributions, wherem, andŝ are the estimated parameters and corresponding standard errors listed in Table 2. Lapatinib concentration absorption from the serum into the tumor was sampled from a Beta(2.3, 1.5) distribution, which has mean 0.61, the average tumor/serum concentration reported in [32].
In silico trials that additionally varied cell motility sampled the diffusion parameter from a Gamma(1, 1) distribution that was truncated below 0.55 and above 4.23, the diffusion parameter thresholds reported in Murray, 2012 [36]. Diffusion parameter values were then divided by 30 to convert from mm 2 /month to mm 2 /day.

In vitro tumor cell pharmacodynamics
We first modeled the in vitro cell pharmacodynamics during lapatinib treatment using a logistic ordinary differential equation (ODE) model (Eq 1). The choice of a logistic model over simpler exponential models was made based on exploratory analysis of viable cell dynamics over time (Fig 2A and 2B). During treatment with low concentrations of lapatinib, the growth rate of cells leveled off after 5 days; this change in the growth rate as a function of the cell number is better captured by a logistic than an exponential model. Note that the birth, death, and clearance rates and carrying capacity are restricted to be non-negative real numbers. The carrying capacity, K, is assumed to be invariant to lapatinib concentrations and is estimated as the maximum observed number of viable cells on day 5 in the absence of lapatinib. The birth, death, and clearance rate parameters b, d, and c, respectively, are functions of the lapatinib concentration; however, the exact functional form describing the relationship between these rates and the lapatinib dose is unknown.
To estimate the birth, death, and clearance rates that best capture cell growth patterns, we implemented a grid search algorithm (Eq 2) to minimize the squared distance between observed and predicted cell numbers for each experimental setting in order to determine the parameters that best reproduce observed cell numbers. The patterns of birth, death, and clearance rates against lapatinib concentration are shown in S1A-S1C Fig. We observed that, when allowing the clearance rates to vary for different concentrations of lapatinib, only the death rates show a discernible pattern with increasing lapatinib dose (S1B Fig). However, if the clearance rate is constrained to a constant, we observed that both birth and death rates demonstrated clear dose-dependent correlations (S1D-S1F Fig). We thus adopted the latter approach, leading to birth rates exponentially decreasing and death rates exponentially increasing with escalating lapatinib concentrations. To further investigate the effects of a constant clearance rate, we studied a wide range of values for the clearance rate (0.0, 0.1, 0.2, 0.3, 0.4, and 0.5), and found only minor differences in the birth and death rates estimated for different clearance rates. We therefore selected the simplest assumption with a clearance rate of 0, which reduces the final model to Eq (1), with c = 0 in Eq (1b).
We then modeled the relationship between birth and death rates and lapatinib dose as exponential distributions (Eq 3). We estimated the distribution coefficients for birth and death rates using nonlinear least squares regression, which are shown in Table 2. We found that birth and death rates are both significantly associated with lapatinib dose (p birth = 1.204 × 10 −5 , p death = 3.351 × 10 −6 ).
Once we established the relationship between the birth and death rates and the lapatinib concentration, we modeled the growth trajectory for different lapatinib concentrations over time (Fig 2C and 2D) using the model described in Eq 1. We found that the number of viable cells predicted by our model agrees with the observed numbers of viable cells on day 5 ( Fig  2E). The model also predicts a clear negative correlation between the lapatinib concentration and the number of viable cells at day 5; however, our model agrees less well with the observed number of dead cells on day 5 ( Fig 2F). One possible reason for the discrepancy is that the number of observed dead cells in much smaller than the number of observed viable cells, and hence birth and death rates were determined predominantly using the observed viable cells. Note that our approach does not capture the peak in the number of dead cells on day 3; we believe that this peak does not represent true underlying biology but rather represents a technical artifact since, if we were to include the peak in our analysis, we would estimate growth and death rates that are inconsistent with the know action of anti-cancer agents (i.e. increasing doses lead to lower proliferation and/or enhanced death). Designing a model to capture this peak would likely lead to overfitting and hence decreased predictive abilities. For these reasons we decided to use the approach outlined above to describe the cell growth trajectories in the presence of lapatinib treatment.

Patient-derived toxicity constraints
To identify optimum treatment schedules for the delivery of lapatinib to patients, our approach needs to take into account clinically determined toxicity constraints-the maximally tolerated dose (MTD) for a particular time interval that does not lead to dose-limiting side effects. The maximum allowable oral dose as a function of the number of treated days was constructed from two clinical studies: Thiessen et al. [16] and Chien et al. [17]; a third point was selected based on clinical expertise in order to estimate the shape of this function (Fig 3A). The three toxicity limits we then fitted to both a linear (Eq 4) and an exponential (Eq 5) model ( Table 2). We found that the exponential model better explains the observed limiting toxicity as a function of the number of treatment days (AIC values: linear 49.48 vs. exponential 5.712). Limiting doses were extrapolated based on the exponential model for 1 to 5 treatment days in a 5-day treatment cycle.

In vivo lapatinib pharmacokinetics
In addition to modeling tumor cell dynamics in the presence of lapatinib treatment, our approach also requires a model of lapatinib pharmacokinetics-the dynamics of lapatinib absorption into the blood stream, absorption into the tumor, and clearance from the tumor. We first modeled lapatinib absorption into the blood stream using pharmacokinetic parameters from clinical trials ( Table 4 in Burris et al [37] and Table 4B in Chien et al. [17]). We then fit a logarithmic function (Eq 6) to this data to describe the relationship between the limiting oral doses identified above to the maximum plasma concentration (Fig 3B, Table 2).
To convert the plasma concentration to tumor concentration, we used patient pharmacokinetic data reported in Table S4 in Vivanco et al. [32]. In this study, the authors found the mean tumor/serum concentration to be 0.61 in biopsied GBMs following 1 week of daily lapatinib treatment. The resulting serum to tumor conversion function is described in Eq 7.
Finally, we constructed a time-dependent pharmacokinetic model using parameters C max = 4 hours, i.e. the time to reach the maximum serum concentration after drug administration, and T 1/2 = 24 hours, i.e. the time to reach half of the maximum concentration after C max [38]. A linear function was used to model lapatinib uptake from 0 to 4 hours, and an exponential function was used to model the subsequent decay pattern, with a half-life of 24 hours (Fig 3C).

Optimum dosing strategies for lapatinib treatment
Combining tumor cell dynamics with patient derived lapatinib pharmacokinetics provides the in vivo tumor growth trajectory-the basis for comparing different dosing strategies. In particular, we compared five treatment schedules that provide the MTD per day according to the identified toxicity constraints, including one continuous dosing strategy and four pulsatile dosing strategies (Table 3).
Optimum dosing strategies using a logistic ODE growth model. We first investigated the five individual treatment schedules using a logistic ODE growth model (Eq 1), which is the simplest growth model possible under our assumptions. The in vitro viable cell trajectories for the five treatment schedules for both short-term response (one 5-day cycle) and long-term response (twenty 5-day cycles) are shown in Fig 4A and 4B. Based on this analysis we observed that none of the five schedules is capable of reducing the total number of viable cells either in the short-term or long-term. In the short-term, all pulsed dosing schedules with large up-front loading doses result in slower viable cell growth than the continuous dose schedule on the initial high dose days. However, once switching to the treatment holidays, viable cell growth in the pulsed schedules outpaces that of the continuous schedule, resulting in a sharp increase in the numbers of viable cells. In the long-term, pulsed schedules result in large fluctuations in the number of viable cells. In contrast, the number of viable cells in the continuous schedule increases smoothly over time. In addition, the number of viable cells in the continuous dosing schedule converges to a limit that is lower than the carrying capacity whereas the maximum numbers of viable cells in the four pulsatile dosing schedules are higher than that in the continuous dosing schedule. Finally, comparing the number of viable cells at the end of treatment to a control schedule (no treatment), we see that the continuous dosing schedule results in a 60% reduction in the total number of viable cells at the end of treatment compared to no treatment, whereas the best pulsatile treatment has less than a 30% reduction (Fig 4C). Although lapatinib fails to arrest further tumor expansion, our data suggest that lapatinib has the ability to reduce overall tumor burden below the carrying capacity under a continuous dosing schedule (schedule 1). Additionally, continuous dosing of lapatinib decreases the number of viable cells at the end of 20 treatment cycles compared to the pulsatile treatment schedules.
Optimum dosing strategies using a diffusion PDE growth model. Human GBM is known to be highly diffusive and infiltrative [33,36,[39][40][41][42][43][44][45]. Hence our logistic ODE growth model based on in vitro growth may not accurately reproduce some of the important characteristics of human disease as it does not account for tumor cell diffusivity and heterogeneity. In order to better scale our model to human disease, we performed additional analyses using a logistic growth diffusion model (Eq 8). Diffusion models have been widely used to model tumor growth [33,36,[46][47][48][49]; for example, in Swanson et al. [33] the authors showed that they can accurately predict patient survival using a diffusion model for tumor growth, supporting the applicability of such a model for investigating human GBM growth kinetics.
To incorporate the diffusion model into our approach, we used a proliferation rate of ρ = 0.387/day and diffusion parameter D = 0.03mm 2 /day, which represent the patient-derived median proliferation and diffusion rates reported in Murray, 2012 [36]. We let the initial solid and isolated cell tumor radii be 1.5cm and 2.5cm, respectively, where the solid and isolated cell tumor volumes represent tumor cell concentrations of 80% and 16%, respectively, of the maximum cell concentration; 1.5cm was the initial solid tumor radius used in [33], and 2.5cm is in the isolated cell radius range reported in Swanson et al. [33]. Note that we chose 2.5cm rather than the median reported isolated cell radius of 2.84cm to allow simulated tumors sufficient time to grow before reaching a radius of 3cm, which is consistently considered to be the "fatal" radius [33]. Finally, we let the lower threshold for detection of cancer cells on an MRI be 8,000 cells/mm 3 , and the carrying capacity, κ, be 10 6 cells/mm 3 . To model GBM proliferation during treatment, we scaled the proliferation rate by the birth and death rates during treatment estimated from the in vitro data (Eq 9). The resulting diffusion model in the presence of treatment is described in Eq 10.
The long-term tumor growth trajectories for the five MTD treatment schedules are shown in Fig 4D. Consistent with our analysis based on the ODE model, we observed that even when using the PDE model, none of the treatment schedules are able to reduce the total tumor volume. However, as before, the continuous dosing schedule again results in slower tumor growth and a reduced tumor volume at the end of treatment compared to the pulsatile dosing strategies ( Fig 4F).
Next, we simulated a clinical trial testing the five MTD treatment schedules in a cohort of 50 patients. We varied the GBM growth rates and pharmacokinetic parameters across patients (see Methods) and considered the time of death of a patient to be when their tumor reached 3cm in diameter [33]. The resulting Kaplan-Meier survival distributions are shown in Fig 4E. We observed that patients treated with the continuous dosing schedule (schedule 1) have the best survival outcomes, while patients treated with the one treatment day/cycle schedule (schedule 5) have the worst survival outcomes, which is consistent with the predicted growth rate trajectories in Fig 4D. Moreover, using a log-rank test, we found that the differences between the continuous dosing schedule and schedules 3, 4, and 5 are significant (p schedule3 = 0.0316, p schedule4 = 7.7 × 10 −16 , p schedule5 < 2 × 10 −16 ). The difference between continuous dosing and schedule 2 (4 treatment days/cycle) is not significant, which is not surprising because the two schedules are very similar.
We observe that the median survival times predicted by our model (80-90 days for schedules 4 and 5) may underestimate survival times observed in the clinic (12-14 months with agressive treatment). However, our analysis suggests that schedules 4 and 5 are only slightly better than no treatment at all. Thus our results are not unreasonable given that median survival for patients receiving no treatment or only palliative antitumor therapy is 8.8 weeks [50], so we are comfortable using this analysis to compare treatment schedules.
By using the MTD on each treatment day for the individual schedules, the total dose administered during a treatment cycle is different for each schedule because administering more frequent treatments of smaller allotments allows administration of more total dose during a treatment cycle compared to less frequent administrations of large doses. The rationale for this approach is to maximize the total dose for each schedule design while still remaining under the MTD. Alternatively, to investigate whether the differences in continuous vs pulsatile dosing are due to differences in the total amount of drug administered, we also investigated tumor growth during five treatment schedules with a fixed total dose ( Table 4). The long-term growth trajectories for the five fixed total dose treatment schedules are shown in Fig 4G, the Kaplan-Meier analyses from the simulated clinical trial using the fixed total dose schedules is shown in Fig 4H, and the comparison of tumor volumes at the end of each treatment are shown in Fig  4I. When implementing this approach, we observed that the ranking of schedules with regard to their treatment efficacies is consistent with our previous results using the maximum tolerated administered dose per day. However, the fixed total dose schedules result in larger tumors at the end of treatment compared to the MTD schedules. Therefore we performed the remainder of our analysis using the MTD schedules listed in Table 3 as our goal is to identify the best overall treatment strategy regardless of the total drug dose administered in one treatment cycle.
The long-term growth trajectories of tumors with varying pre-existing resistance are shown in S2 Fig, and a comparison of tumor volumes at the end of treatment between the five schedules across different levels of pre-existing resistance is shown in Fig 5A. We found that for 0% and 1% pre-existing resistance, continuous dosing leads to better efficacy than any of the pulsatile dosing schedules; however, when there is 5% or 10% pre-existing resistance, the continuous dosing and 4 doses/cycle schedule (schedule 2) perform similarly, though still better than the other pulsatile dosing schedules. These results indicate that the continuous dosing schedule is preferable to pulsatile dosing even in the presence of a resistant cell population, although the difference among the individual treatment schedules becomes smaller at higher levels of preexisting resistance.
Non-uniform treatment delivery to the brain. A major reason for GBM treatment failure is insufficient drug delivery to the brain. The blood brain barrier prevents most systematically-delivered molecules from being transported into the brain [67]. Indeed, the tumor concentrations of lapatinib in the patient cohort from Vivanco et al. [32] vary considerably, ranging from 70-3826 nM, as do the plasma/tumor ratios, and the blood brain barrier may be one reason causing such variable lapatinib concentration reaching the tumor.
To establish if the blood brain barrier affects lapatinib treatment response, we made the assumption that the blood brain barrier causes less drug to enter the tumor leading to a decreased drug effect. We adapted our diffusion model to account for reduced drug delivery to the brain by introducing a factor Q in the proliferation rate function (Eq 13). We use this model to assess how much a 10%, 25%, and 50% decrease in the lapatinib concentration reaching the tumor from the serum affects our simulated treatment outcomes. The growth trajectories resulting from different amounts of lapatinib reaching the tumor for each treatment schedule are shown in S3 Fig. We also compared tumor volumes at the end of treatment between treatment schedules across different levels of lapatinib reaching the tumor ( Fig 5C). As expected, we observed that decreasing the amount of lapatinib that crosses the blood brain barrier decreases the effect of treatment. However, the decreased effect of treatment is similar across all treatment schedules; hence continuous dosing still remains the best treatment option for reducing longterm tumor growth, even when there is reduced availability of lapatinib in the brain. The effects of motility on lapatinib treatment response. Tumor cell motility has been described as a source of variability in treatment responses among a patient cohort and measures of motility can be predictive of treatment response [36]. In order to asses the effects of tumor cell motility on treatment outcomes, we simulated tumor growth using Eq (8) for varying values of the diffusion parameter, D, estimated from patient data (ref). For our analyses we utilized D equal to 0.01833, 0.0333, 0.0667, 1.333, and 1.667 mm 2 /day, where 0.01833 mm 2 / day is the lower limit and 1.667 is above the upper limit reported in Murray, 2012 [36]. Long-  Fig  5C). Although we found, as expected, that tumors with higher motility grow faster during all treatment schedules, the continuous dosing schedule performs better than all pulsatile dosing strategies across all diffusion parameters tested. We then investigated whether the relationship between motility and tumor volume (in mm 3 ) at the end of treatment depends on the treatment schedule using linear regression with an interaction term between the diffusion parameter and the schedule (Eq 17). We found that schedules 3, 4, and 5 are significantly more affected by increased motility than schedule 1 ( Table 5), indicating that the effectiveness of continuous dosing is less comromised by increased cell motility than that of a pulsatile schedule with only 1, 2, or 3 treatment days per cycle. This observation suggests that a continuous dosing strategy is more effective than pulsatile dosing even with a large extent of cell motility, consistent with our previous results.
Tumor motility heavily depends on the location of the tumor in the brain, as malignant glial cells diffuse more quickly in white matter than gray matter [68], [69]. Thus, tumor motility can be highly variable between patients, depending on where their tumors are located. To further investigate how variability of motility across a patient cohort affects treatment outcomes and survival, we simulated a clinical trial of the five MTD treatment schedules with a varying diffusion parameter in addition to varying cell dynamics and pharmacokinetics parameters (see Methods). The resulting Kaplan Meier analysis is shown in Fig 5D. We observed that patients treated with a continuous dosing schedule have the longest survival, whereas patients treated with schedule 5 have the worst survival outcomes, consistent with our previous results. Using a log-rank test, the difference between the continuous dosing schedule and schedules 4 and 5 are significant (p schedule4 = 0.0104, p schedule5 = 1.33 × 10 −6 ).
Incorporating the Go-or-Grow mechanism into our mathematical framework. Glioblastoma is known to be highly hypoxic [34,70], and tumor cell invasion is one of many different rescue mechanisms that allow tumors to survive hypoxic conditions [34,71]. According to the Go-or-Grow hypothesis, in low oxygen conditions, tumor cells undergo a phenotypic switch from proliferative to invasive/migratory phenotypes [72]. Multiple mathematical models have been developed to describe the effects of this phenomenon on tumor growth [35,[73][74][75][76]. We adapted Model II described in [35] to investigate the five laptinib treatment schedules under the assumption of the Go-or-Grow model. Our model is decribed in Eqs 14-16 and parameterized using experimentally derived values, as in [35]. We used vascular density v ¼ 1 2 ; oxygen concentration in normal tissue σ 0 = 2.068 × 10 −9 mol mm −3 ; proliferative to migratory switching parameter λ 1 = 4.134 nmol mm −3 (twice the normal tissue concentration); migratory to proliferative switching paramter λ 2 = 0.5, 1, and 2; oxygen diffusion rate D σ = 1.51 × 10 2 mm 2 day −1 ; oxygen supply rate h 1 = 3.37 × 10 −1 day −1 ; and oxygen consumption rate h 2 = 10 −2 mm cell −1 day −1 (in the range of values reported in [35]). When implementing this model, we found that continuous dosing performs better than all pulsatile dosing schedules (S5 Fig), consistent with our results above.

Discussion
Here we have presented a computational analysis platform of human GBM growth and treatment response, parameterized using in vitro tumor cell data measured based on the SF268 GBM cell line during lapatinib exposure. Our data demonstrate that lapatinib concentrations are negatively correlated with cell birth rates, but positively correlated with death rates. However, the necessary lapatinib concentration needed to arrest further tumor expansion is calculated to be 2,180nM in vitro, a concentration that cannot be reached in GBM patients based on current toxicity limits and oral absorption profiles. To achieve this concentration, our absorption model predicts that an oral dose of 34,000mg is required-5 times greater than the highest dose ever tested. The therapeutic potential seems to be restricted by the transport of lapatinib from the gastrointestinal tract to tumor sites. Due to the logarithmic relationship between oral intake and serum concentration, a higher than tolerated dose is required to achieve tumor reduction. Other methods of delivery, such as intravenous injection, may need to be considered to bypass the current delivery problem. This result agrees with the conclusions from early clinical trials that lapatinib failed to reduce tumor sizes [15], [16]. Despite these negative results, our model suggests that continuous lapatinib dosing does have positive effects. Continuous dosing was able to reduce long-term tumor burden significantly more than pulsatile treatment schedules. These findings were consistent with incorporating intratumor heterogeneity due to the presence of drug-resistance cells, variable penetration of the blood brain barrier leading to a reduction of lapatinib concentration in the tumor, and tumor cell motility. Optimization of dosing strategies for cancer treatments is a complex problem. Our mathematical model provides a quantitative relationship bridging in vitro tumor dynamics and in vivo lapatinib dose and schedule, taking into account the blood brain barrier, GBM motility, heterogeneity, and phenotypic changes that occur in hypoxic microenvironments. Our approach allows for the systematic comparison between different treatment strategies and their effects on tumor growth. Although none of the tolerable doses and schedules evaluated were found capable of arresting tumor expansion, our approach predicts that continuous treatment within these parameters is likely to be more successful at slowing down tumor growth. Further optimization of effective dosing strategies using our computational platform will benefit from additional dose-limiting toxicity studies that incorporate multiple pulsatile schedules.

S1 Fig. Birth, death and clearance rates for different concentrations of lapatinib. A-C:
There is no obvious association between birth and death rates with lapatinib concentrations when clearance rates are allowed to vary. D-F: Constraining the clearance rates to be zero, birth rates decrease with increasing lapatinib concentrations and death rates increase lapatinib concentrations. The solid black lines show the relationships between birth and death rates and lapatinib concentration assuming an exponential function. The exponential function is selected to ensure that birth and death rates are always positive for all concentrations of lapatinib.