Quantifying glioblastoma drug response dynamics incorporating resistance and blood brain barrier penetrance from experimental data

Many drugs investigated for the treatment of glioblastoma (GBM) have had poor clinical outcomes, as their efficacy is dependent on adequate delivery to sensitive tumor cell populations, which is limited by the blood-brain barrier (BBB). Further complicating evaluation of therapeutic efficacy, tumors can become resistant to anti-cancer drugs, and it can be difficult to gauge the extent to which BBB limitations and resistance each contribute to a drug’s failure. To address this question, we developed a minimal mathematical model to characterize these elements of overall drug response, informed by time-series bioluminescence imaging data from a treated patient-derived xenograft (PDX) experimental model. By fitting this mathematical model to a preliminary dataset in a series of nonlinear regression steps, we estimated parameter values for individual PDX subjects that correspond to the dynamics seen in experimental data. Using these estimates, we performed a parameter sensitivity analysis using Latin hypercube sampling and partial rank correlation coefficients. Results from this analysis combined with simulation results suggest that BBB permeability may play a slightly larger role in therapeutic efficacy than drug resistance. Our model and fitting technique to estimate parameters from data may be a useful tool in aiding further exploration of these challenges in future studies of drug efficacy with larger datasets.

analysis combined with simulation results suggest that BBB permeability may play a slightly larger role in therapeutic efficacy than drug resistance. Our model and fitting technique to estimate parameters from data may be a useful tool in aiding further exploration of these challenges in future studies of drug efficacy with larger datasets.

Introduction
Glioblastoma (GBM) is an aggressive primary brain cancer that is notoriously difficult to treat due to its diffuse infiltration into surrounding normalappearing brain [10]. These diffusely invading GBM cells cannot be completely resected surgically [2], and are difficult to target with radiation therapy while sparing normal brain [5]. As a result, clinicians rely on chemotherapy to treat the full extent of the tumor. However, chemotherapeutic efficacy can be limited in two main ways: there may be insufficient delivery across the blood-brain barrier (BBB), and the tumor may develop resistance to therapy.
The BBB acts to keep pathogens and many toxins out of the sensitive brain tissue. Angiogenesis in dense tumor regions induces disruption of the BBB, potentially allowing chemotherapeutic drugs to "leak" into these tumor regions. Current dogma in neuro-oncology holds this as being largely sufficient to treat the tumor, but GBM cells invade beyond these regions into tissue where the BBB remains rather intact [19]. Further, tissue interstitial pressure and drug properties such as lipophilicity and polarity may influence the delivery of drugs across angiogenesis-induced BBB "leaks" [16]. Due to these factors, it remains unclear whether the delivery of BBB-impermeable antineoplastic agents reaches adequate concentrations throughout the tissue to provide the anticipated therapeutic effect.
Drug insensitivity or resistance is also a key suspect behind unsuccessful results treating GBM with molecularly-targeted therapies [7,20]. GBMs frequently present with gene mutations or amplifications for a number of targets, such as epidermal growth factor receptor (EGFR), for which therapies already exist [3,4,8,15,18]. Due to the spatial heterogeneity of GBM, it is possible that these targets have been identified for a subpopulation of cells predominant in the dense core of the tumor where biopsies were taken, but are in fact less common in the invading portions of the tumor. This would result in a significant population of cells being treated with drugs they are not sensitive to, which could explain why trials with such targeted agents have failed [1,6,11,17]. It is notable, however, that many of these drugs were not developed specifically for brain, but for cancers elsewhere in the body and subsequently tried in brain cancer, again raising the question of adequate delivery across the BBB.
"Specifically for EGFR or other kinase-targeted inhibitors, emergence of compensatory signaling pathways may lead to acquired resistance to ther-apy." In order to explore both the contributions of inadequate delivery of therapy across the BBB and drug resistance or insensitivity, we developed a minimal mathematical model based on preclinical experimental data. First, we describe model development based on this data and steps to estimate parameter regimes via data-fitting. Next, we explore the global model parameter sensitivity to understand how these parameters impact model outcomes. Finally, we run model simulations for the data-derived parameter regimes to assess the relative contributions of drug distribution and sensitivity, and discuss how it might be useful in assessing results from future experiments comparing different PDXs or different drugs. Overall, our model suggests that the influence of drug permeability may be more impactful than the degree of resistance for a given baseline sensitivity. Thus, in order to improve treatment outcomes, it is critical to determine predictors of drug distribution in individual patients' tumors and surrounding brain tissue to ensure invading tumor cells are adequately exposed to the therapy.

Treatment Exposure and Sensitivity Model
Our ordinary differential equation (ODE) model of tumor growth and treatment response accounts for both variable treatment exposure and differential sensitivity to treatment by different tumor subpopulations. Development of this model was informed by experimental observations, which were also used to determine relevant parameter regimes for running simuliations.

Experimental Data
The form of our model was based on experimental data testing an EGFRtargeted antibody drug conjugate (ADC) [14]. These experiments were performed using a patient-derived xenograft (PDX) model of GBM, in which patient-derived GBM cells are implanted in rodents [9,13]. The growth of these preclinical tumors was monitored via bioluminescent imaging (BLI, Figure 1). Since BLI flux is linearly correlated with tumor cell number [12], this provided us with a close approximation of tumor cell populations across time. Importantly, the PDX tumors were implanted both in the flank, where there is no BBB, and in the brain, which allowed us to compare treatment effect in tumors with and without BBB impediments to drug distribution. Moreover, there were control groups in both sites (flank and intracranial) that received no treatment.
In the untreated groups, BLI flux increased exponentially over time, indicating exponential growth in the total tumor cell population. Treated groups showed similar exponential growth in the absence of treatment at early time points, followed by a precipitous decline with the initiation of therapy (begun once tumors reached a set volume). This period of decline Colors represent the BLI radiance in photons per second per area in cm squared.., which is related to the BLI flux (measured in photons per second) for the total area. was relatively short-lived, and after that BLI flux increased exponentially again, albeit at a slower rate. The precipitous decline in experimental tumor volume followed by a trajectory of slower regrowth than that prior to therapy suggested the presence of two distinct tumor cell populations: one that responded to therapy (sensitive, denoted s), and one that continued to grow in spite of therapy (resistant, denoted r). These cell populations are assumed to proliferate at the same rate, ρ, differing only in their responses to treatment. The ADC therapy (denoted A) was delivered every seven days at dose A dose , and decays at rate λ. In response to therapy, which achieves a fraction γ of the dosed concentration in the tumor tissue, cell death is induced in each population at rates (µ s,r ), respectively. These dynamics are schematized in Figure 2.

Model Equations
Our model consists of three coupled ordinary differential equations describing the dynamics of both cell populations (s, r) and the ADC (A): where parameters and their definitions are outlined in Table 1, and their derivations can be found in Section 2.3.
In the absence of the ADC, both sensitive and resistant tumor populations grow exponentially, at proliferation rate ρ. However, the two populations differ in sensitivity to the ADC, A, which is captured by the druginduced apoptosis rates µ s and µ r (for sensitive and resistant populations, respectively). The terms for tumor cell death due to ADC are further modified by factor γ, which represents the proportion of cells exposed to ADC. We assume that the ADC is readily distributed to flank PDXs such that tumor cell exposure is high (γ = 1), but that the BBB limits this distribution for intracranial PDXs (0 ≤ γ ≤ 1). In order to capture the ADC dynamics, we let A dose (n) represent the nth dose, with doses administered every seven days, as noted by the dirac delta function δ(t − 7n). The ADC then decays at rate λ. Table 1: Model Parameter Definitions and Values. Parameter value ranges were estimated through fitting the model to experimental data or parameters were confined to a value range by their theoretical meaning, except in the case of ADC parameters, as described in Section 2.3.

Symbol Definition
Value Range Units ρ cellular proliferation rate 0.2 to 0.5 day −1 µ s ADC-mediated sensitive cell kill rate 1 to 10 proportion of tumor exposed 0 to 1 unitless A dose ADC given in a single dose 0.1 mg This model can be solved analytically, as shown in Appendix A. For simplification, at any given time t, C(t) represents the total number of cells, calculated by the sum of sensitive s(t) and resistant r(t) cells. This total cell number was used in Section 2.3 for comparing with bioluminescence imaging data, which shows the total tumor cell population. The initial resistant proportion of total implanted cells is denoted by q = r 0 /C 0 . Similarly, the extent to which resistant cells r are less sensitive to the agent than the sensitive cells s is denoted by the ratio z = µ r /µ s , which is bounded between 0 and 1 to ensure that µ r is a fraction of µ s in the regression-based parameterization in Section 2.3. With these notational changes, we can then write the analytical solution (derived in Appendix A) as where Utilizing this solution (2), the model can be parameterized through comparison of simulations to temporal BLI data.

Data-Based Parameter Estimation
Most model parameters were unknown, with the exception of ADC-specific parameters: the timing of dose administration and dose amounts (A dose (n)), as well as the half-life of the drug, which allowed us to solve for the drug decay rate (λ). Dose amounts were adjusted for the weight of each animal (5 mg/kg), so we applied the average initial animal weight of 20 g to obtain the constant ADC dose, A dose = 0.1 mg used in simulations. All of the remaining model parameters were determined through several iterations of fitting the model via least squares regression to preliminary BLI data from an experiment. The various arms of the experiment included untreated and treated groups of subjects, as well as flank and intracranial tumor sites to separate out BBB influences. By fitting the model to these various subgroups, we are able to identify and estimate each of the parameters, as described below.
Step 1: Fit to untreated data to estimate growth rate, ρ. When fitting the model to untreated data, since the ADC is not injected (A = 0), the model's treatment components zero out and only an exponential growth function remains: C(t) = C 0 e ρt . Fitting this model function to untreated data via least squares regression with the lsqcurvefit function in Matlab R (MATLAB Release 2018b, The MathWorks, Inc., Natick, Massachusetts, United States) (Figure 6), we were able to obtain estimates of the tumor proliferation rate ρ and the number of viable implanted PDX cells C 0 . (Note that while a consistent number of cells are initially implanted for each subject, C 0 is in fact unique for each, as a variable number of cells die off, possibly due to an inability to establish themselves in the proper microenvironment for growth.) This yielded subject-specific values for ρ and C 0 , and the mean ρ was recorded as the net proliferation rate for the cells of the particular PDX line used in the experiments grown in either the flank or intracranial setting.
Step 2: Fit to treated flank data to estimate µ s , z, and q. Using the estimated net proliferation rate, ρ, from the previous step, we proceed to fit the treated data in the flank. We assume the flank-specific estimate of ρ remains the same in the treated case as untreated, since the microenvironment remains similar and any effects on growth can be encapsulated in the treatment effect term. Additionally, since the tumor was injected in the flank, there is no BBB effect to limit the proportion of tumor exposed to the ADC, such that the exposure parameter γ = 1. Pairing this with other known parameters (see Table 1), the only remaining unknown parameters are the cell death rates due to drug, µ s and µ r , and proportion of implanted cells that are resistant, q. Using the definition z = µ r /µ s and the analytical solution of the model (2), we can then apply a nonlinear least squares regression (again using lsqcurvefit) to fit subject-specific parameters for parameters µ s , z, and q ( Figure 6).
Step 3: Fit to treated intracranial data to estimate γ, z, and q. Proceeding to fit the data from treated intracranial tumors, we apply the same approach to estimate parameters as in the flank, this time assuming that the intracranial-specific estimate of ρ from the untreated setting remains the same for the treated intracranial tumors due to a similar microenivronment. We further assume that the cell death rate due to therapy for the sensitive tumor subpopulation, µ s , is the same intracranially as in the flank setting and estimate parameter γ, the fraction of tumor exposed to therapy. At the conclusion of these steps (summarized in Figure 3), all unknown model parameters had net and individual estimates (listed in Table 1). Using these values then allowed us to run simulations in a reasonable range of parameter values, as well as to perform a model sensitivity analysis to understand how variability in these values affect model outcomes.

Parameter Sensitivity Analysis
Due to the uncertainty and variability in our parameter estimates, it was important to better characterize the effects of parameters on model results.
To do this, we conducted a parameter sensitivity analysis via Latin hypercube sampling (LHS) and partial-rank correlation coefficients (PRCC). To perform the LHS analysis, we first drew 1000 equiprobable samples for each unknown parameter, including the initial condition C 0 , from a statistical distribution of values. These distributions were informed by our fits of the preliminary data when available; in the case of the unitless parameters, we assumed a uniform distribution on the interval [0, 1]. These samples were then randomly paired in a Latin hypercube scheme to run a series of 1000 Monte Carlo simulations. Using these simulation results, we then computed PRCCs between each parameter and two different model outcomes across all time points: the total number of tumor cells and the fraction of tumor that is resistant (Figure 4). Further details for about this method and the code files used are available on GitHub: https://github.com/scmassey/modelsensitivity-analysis.
3.1.1 Total tumor population depends most strongly on proliferation rate, followed by treatment response parameters.
At early time points, particularly before the initiation of therapy, the tumor population is strongly positively correlated with both the initial number of cells implanted, C 0 and proliferation rate ρ ( Figure 4A). However by 30 days, or after approximately three doses of therapy, the population is strongly positively correlated with ρ and the effect of C 0 begins to wane. At the same time, drug sensitivity of the s cell population, µ s , and exposure to drug, γ are strongly negatively correlated. Resistance factor z, which determines the fraction of drug sensitivity remaining in the r cell population, is also negatively correlated, but less strongly, and only reaches a PRCC value of −0.5 after 100 days. This suggests that resistance z is less impactful than the baseline sensitivity of the tumor population and the distribution of drug in the tissue. As expected from our substitution of µ r = zµ s , which results in the coefficient −γzµ s in the term describing drug induced apoptosis for the equation describing the r population (1b), which mirrors that for the s population (1a), −γµ s , parameters γ and µ s track together in the sensitivity analysis (overlapping lines in Figure 4). Thus, sensitivity analysis is unable to compare the differential impacts of these two parameters. However, we note that since we had preliminary data in both the treated flank as well as the treated intracranial PDX settings, we were able to obtain parameter estimates for these by keeping γ = 1 in the flank setting, and assuming that µ s is the same intracranially as in flank.

Figure 4: Partial rank correlation coefficients (PRCC) of parameters with respect to (A) tumor cells and (B) fraction of resistant cells (r/(s + r)
= r/C), visualized across simulation time.

Resistant fraction of tumor driven by initial resistant proportion, followed by treatment response parameters.
Prior to the initiation of therapy, only parameter q, the fraction of initially implanted cells that are resistant, is correlated with the proportion of total tumor that is resistant ( Figure 4B). Once treatment is initialized, q remains highly positively correlated, and this correlation decreases slightly over time during the course of treatment. Three other parameters show correlation with the fraction of tumor that is resistant following initiation of therapy, all of which involve drug response. Parameters γ and µ s , representing the degree of exposure to drug and baseline sensitivity of the cells, respectively, are both positively correlated and track together, while parameter z, representing the fraction of treatment sensitivity remaining in resistant cells, is negatively correlated with the fraction of tumor that is resistant. Further, the PRCC values do not vary over the time of the simulation after treatment is initiated and sustained. These correlations are consistent with expectations from the behavior of the system described by the model.

Simulation Results
To more fully explore the effect of parameters on model predicted outcomes, we ran simulations for varied values of the parameters relating to treatment response: γ, µ s , and z, or exposure to therapy, sensitivity to therapy, and degree of reduced sensitivity to therapy in resistant cells, respectively. Codes used to run simulations and plot the results may be found on GitHub: https://github.com/scmassey/treatment-exposure-sensitivity-model.

Treatment exposure is especially important for tumors with lower baseline sensitivity
Comparing simulation results across a range of values for parameters γ and z, we see that γ plays a larger influence on total tumor cells for a given baseline sensitivity than does z. That is, our model suggests that exposure to therapy is more impactful on total tumor burden than the degree of resistance, which is consistent with the parameter sensitivity results of Section 3.1.1.
The simulations also indicate that for higher sensitivity, treatment can be effective at lower levels of exposure ( Figure 5, compare heatmaps for µ s = 3 vs µ s = 7).

Same total tumor burden, different fractions of sensitive vs resistant subpopulations
Further, we observed that there can be distinct differences in the dynamics of the individual cell populations underlying simulations that show the same resulting overall tumor cell population level ( Figure 5A,B). Looking at long time scales (> 50 days)-in this case at 12 weeks or 84 days, the average survival time of the treated subjects-we are able to observe the effect of an extended time of treatment in the simulations. It is particularly notable that one simulation retains a large proportion of treatment sensitive cells ( Figure 5A), while another is made up almost entirely of treatment resistant cells ( Figure 5B). These two simulations share the same baseline sensitivity, µ s , and the same level of exposure, γ, differing only in the degree z to which the resistant subpopulation is insensitive to the ADC.

Parameter Estimation for comparing across subjects and experimental groups
While only performed for a small set of data, we note that our parameter estimation procedure may be useful for comparing results within and across experiments. In this small data set, the variation in growth curves between the five treated flank and intracranial subjects was associated with differences in estimated parameter values. (This data, with overlaid simulations using the fitted parameter estimates are shown in Figure 6.) Comparing

Discussion
Our Treatment Exposure and Sensitivity model describes tumor growth and treatment response incorporating the effects of the BBB on tumor exposure to therapy as well as differential therapeutic sensitivity. It is important to note that tumor heterogeneity may actually provide for many subpopulations with varying levels sensitivity to therapy; we assumed that these cluster towards more or less sensitive, reducing these to two subpopulations. Additionally, the therapy may have an effect not only through induction of tumor cell apoptosis, but could also reduce tumor proliferation. Because our data does not provide for an ability to distinguish between these effects, both are encapsulated in the µ s,r therapeutic effect parameters. This assumption allows us to use the estimates of ρ obtained through fitting the model to untreated flank and intracranial experimental groups in the treated flank and intracranial settings. Parameter sensitivity analysis of the model highlighted the necessity of having the flank and intracranial treatment groups for practical identifiability in obtaining estimates for γ and µ s . This "tradeoff" between therapeutic sensitivity (µ s ) and exposure (γ) was also observed in simulations and is expected given our model formulation. However, the emergence of this dynamic in the creation and parameter estimation of our model underscores that this relationship should be carefully considered in the design of experimental studies for new glioma therapies.
Sensitivity analysis also revealed that parameter ρ is the strongest positive influence on total cell population, as expected, and distinguished between the effect of reduced therapeutic sensitivity (z) and exposure to therapy (γ) in reducing the total tumor population. Not only is there a difference in the magnitude of correlation between parameters γ and z with total tumor burden, there is also a difference in the temporal dynamics of the change in these correlations over time. Further, most correlation coefficients between parameters and total tumor burden change most dramatically at earlier time points, and after approximately 50 days, change relatively little. This suggests that experiments conducted to examine the relationship between exposure and sensitivity to therapy should focus on collecting time course data more densely for the first 7 weeks as compared to longer times. It is also notable that the correlation coefficients between parameters and the resistant fraction of the tumor is quite stable over time. Simulations confirmed the greater impact of therapeutic exposure (γ) relative to reduction in sensitivity (z) foreshadowed by the sensitivity analysis. More importantly, they revealed that for a particular tumor burden at any single time, there can be several parameterizations that fit the data but correspond to different fractions of resistant cells. This indicates that time course data is essential for detecting differences in the contributions of γ and z.
Finally, because the model we have presented is minimal and reduces mechanisms down to a few key parameters, it has greater utility for fitting experimental data to estimate these parameters for individuals. Having these individual parameterizations is key to understanding the extent to which drug exposure and resistance each contributed to variations in outcome. We anticipate using this approach in the future for comparing experimental groups with additional PDX lines as well as other therapies, characterizing the differences in outcome results between those groups. This will enable us to determine whether there is any consistency in parameter values within groups, which may lead to improved understanding of patterns behind drug failures and identification of tumor characteristics that might suggest candidates who would benefit most from an emerging therapy.

A Analytical Solution to Model
Here we present the derivation of the analytical solution for model system (1). Recall that the system consists of the following differential equations, which describe the growth of two tumor cell populations with differential drug sensitivity, s and r, and drug level, A: whereÂ is a vector containing the consecutive drug doses administered, such thatÂ(n) is the dose on the nth pulse (of which there are at most N drug pulses), and 7n denotes the time of the nth pulse (every 7 days). The parameters are given in Table 1 of the main text.
Because the equation for A in the system (4) is independent of s and r, we first solve (4c) for A. First, observe the following rearrangement: Notice that we can use the method of integrating factors, introducing a factor of e λt : so that Integrating the left hand side, and moving the integral inside the sum on the right hand side: Then, integrating the right hand side, we have where θ is the heaviside function.
Now, since λ = ln(2)/7 (i.e., the half life is 7 days), we have the coefficient e 7nλ = e 7n ln(2)/7 = e ln(2)n = 2 n , so the right hand side of (11) becomes Rearranging, we have our solution for A: Now we solve for s and r, and use the fact that (4a) and (4b) differ only by parameter µ s versus µ r to obtain both simultaneously: Integrating both sides: ln(s) = ρt − γµ s A(t)dt + c 1 , where c 1 is a constant of integration to be determined from the initial condition.
Since the solution for s in (18) depends on A(t)dt, we must integrate (14): To use integration by parts, such that u v = uv − uv , choose Then we have: where c n is a constant of integration. Thus, (20) becomes: where c 2 is the sum of all the c n .

20
The solution for r is similar: Now, we assume that the s cells are more sensitive to drug than the r cells, such that µ s > µ r . Defining z to be the ratio z = µ r /µ s , we can substitute µ r = zµ s , and in our fitting process, solve for this degree of differential sensitivity between the s and r cell populations.
Additionally, for fitting purposes, we can combine the cells into one total population: C = s + r. To do this, we write the initial condition C 0 = s 0 + r 0 , and let s 0 = qA 0 such that q represents the proportion of C 0 that is made up of cells that are more drug-sensitive. This gives: Because this is a lengthy expression, in the main text we write C(t) = C 0 e ρt qe −γzµs A(t)dt + (1 − q)e −γµs A(t)dt , where Adt denotes the summation in (20).