Complementary role of mathematical modeling in preclinical glioblastoma: differentiating poor drug delivery from drug insensitivity

Glioblastoma is the most malignant primary brain tumor with significant heterogeneity and a limited number of effective therapeutic options. Many investigational targeted therapies have failed in clinical trials, but it remains unclear if this results from insensitivity to therapy or poor drug delivery across the blood-brain barrier. Using well-established EGFR-amplified patient-derived xenograft (PDX) cell lines, we investigated this question using an EGFR-directed therapy. With only bioluminescence imaging, we used a mathematical model to quantify the heterogeneous treatment response across the three PDX lines (GBM6, GBM12, GBM39). Our model estimated the primary cause of intracranial treatment response for each of the lines, and these findings were validated with parallel experimental efforts. This mathematical modeling approach can be used as a useful complementary tool that can be widely applied to many more PDX lines. This has the potential to further inform experimental efforts and reduce the cost and time necessary to make experimental conclusions. Author summary Glioblastoma is a deadly brain cancer that is difficult to treat. New therapies often fail to surpass the current standard of care during clinical trials. This can be attributed to both the vast heterogeneity of the disease and the blood-brain barrier, which may or may not be disrupted in various regions of tumors. Thus, while some cancer cells may develop insensitivity in the presence of a drug due to heterogeneity, other tumor areas are simply not exposed to the drug. Being able to understand to what extent each of these is driving clinical trial results in individuals may be key to advancing novel therapies. To address this challenge, we used mathematical modeling to study the differences between three patient-derived tumors in mice. With our unique approach, we identified the reason for treatment failure in each patient tumor. These results were validated through rigorous and time-consuming experiments, but our mathematical modeling approach allows for a cheaper, quicker, and widely applicable way to come to similar conclusions.

Glioblastoma (GBM) is the most common primary brain malignancy and is aggressive, 2 heterogeneous, and diffusely-invasive. Despite maximal surgical resection, radiation 3 therapy, and chemotherapy, patients with GBM have a two-year survival of 4 approximately 25%, and five-year survival under 10% [1][2][3]. In efforts to further 5 improve overall survival of patients with GBM, many targeted therapies and other 6 chemotherapeutic drugs have been developed and appear promising in preclinical 7 trials [4]. Yet, many of these drugs fail at the clinical level. To this day, none have 8 surpassed the standard of care established over a decade ago, leaving the median overall 9 survival at a dismal 14.6 months [2]. 10 Drug delivery to the brain is an inherent challenge posed by the blood-brain barrier 11 (BBB), which largely restricts hydrophilic molecules and large macromolecules in the 12 bloodstream from crossing into the brain [5]. While this is a critical aspect of protecting 13 this vital organ against various infections, it also prevents a vast majority of oral or 14 intravenous drugs from being delivered to brain tumors. This inherent limitation in  However, as it is not easy to determine where drug was actually delivered in a human 25 brain, the gold standard for investigating drug failure is via in vivo experimentation 26 with cell lines and animal models. To form definitive conclusions, numerous repeats are 27 required, with tumors grown outside of and within the brain, and the resulting tissue 28 undergoing DNA, RNA, and protein sequence analyses. Recently, Marin et al. did just 29 this for depatuxizumab mafodotin (Depatux-M, ABT-414), an EGFR-directed antibody 30 drug conjugate (ADC) [6]. In brief, they performed intracranial experiments for seven 31 different cell lines, and through numerous experiments with each line, they determined 32 that only two cell lines were responsive to the therapy. Of the remaining five cell lines, 33 two acquired resistance to the drug and three were not responsive to the drug due to a 34 relatively intact BBB. This work was the first of its kind to go so deep and to fully to treatment [6,7]. This is significant because it enables researchers to quickly assess 40 cell lines and determine which deserve deeper experimentation. By reducing time and 41 costs, we hope that this could focus drug development research and increase the speed 42 of acquiring clinically meaningful answers about drug effectiveness.

43
In this paper, we demonstrate that our model, using relatively minimal experimental 44 data, arrives at similar conclusions to the highly sophisticated experiments conducted in 45 Marin et al. for a subset of the cell lines [6]. First, we describe the mathematical model 46 as presented in Massey et al. [7]. Then, we provide a brief description of the cell lines 47 and experiments that were used. We briefly discuss the methodology for Xenograft (PDX) National Resource [8] for their known differential response to 54 Depatux-M (see S1 Table for further details). These cell lines were transduced with 55 firefly luciferase (F-luc) to allow for bioluminesence imaging (BLI). As the resulting BLI 56 flux is linearly proportional to the total number of luminescing cells, we can 57 non-invasively monitor tumor progression in vivo over time [9]. In brief, patient-derived 58 GBM cells were implanted heterotopically into ten athymic nude mice per cell line.

59
After the tumors were established, PDXs were assigned to therapy groups, including a In this paper, we utilize the Treatment Exposure and Sensitivity model, first presented 84 in Massey et al. [7]. It was designed to differentiate the relative contributions of cell 85 sensitivity and drug delivery to the overall tumor response to therapy. We used the 86 previously described preclinical experimental data to parameterize this model.

88
The model is a system of three ordinary differential equations (1) capturing the 89 dynamics of the drug (A) and two underlying tumor populations: one subpopulation 90 that is highly sensitive to the treatment (H) and another that is less sensitive (L). The 91 model, as listed below, has the analytical solution Parameters for the model are summarized in Table 1. were calibrated in a sequential process enabled by the experimental design [7].

101
The calibration of the proliferation rate ρ was done in two parts for both the flank and 102 intracranial settings. First, as the sham control experiments have no therapeutic effect, 103 the model can be evaluated without treatment, thereby simplifying to an exponential Natick, Massachusetts, United States) lsqcurvefit function, the logged-data can be 106 used to fit the initial number of viable tumor cells and the proliferation rate ρ. Since 107 proliferation rate is both intrinsic to the PDX line and specific to the microenvironment, 108 the subject-specific fits obtained in fitting the sham control data was used to develop an 109 informative Gaussian prior. This was then used to estimate the proliferation rate for 110 treatment in the same microenvironment.

113
The calibration of the sensitive and insensitive proportion q and respective rates µ H , µ L 114 was specific to the flank or intracranial settings. Since the flank does not have the 115 limitation of the BBB, we can assume that the entirety of the tumor is exposed to the 116 therapy and that γ = 1. The remaining parameters q, µ H , and µ L were sampled from 117 uninformative priors spanning the ranges in Table 1. The selected parameter set was 118 subject-specific and minimized the sum-of-squared-error for that subject in the flank.

119
With five subjects in the flank, we formulated a Gaussian prior for the µ H sensitivity 120 using the selected parameter sets. Using this prior, we proceed with a similar approach 121 to select remaining parameter sets in the intracranial setting. 122

123
Applying this fitting algorithm to each of the three PDX lines, variations in the model 124 parameters captured the wide range of treatment response dynamics in the data, both 125 across and within experiments. The resulting fits are displayed in Fig 2,     In contrast, estimates of the less-sensitive cell kill rate (µ L ) captured the differences 142 between PDX lines and was consistent across subjects. In the flank experiments, the 143 average rates for GBM6, GBM12, and GBM39 were 2.2, 1.7, and 0.6 mg −1 day −1 , 144 respectively ( Fig 3A). This suggests that, of the PDX lines investigated, the  Fig 3B). This suggests that the less-sensitive population L of GBM39 appears to be 150 less sensitive than the GBM6 and GBM12 intracranially.

151
Proportion of tumor exposed to drug (γ) varies by PDX line

152
Estimates of the proportion of tumor exposed to the drug (γ) also captured the 153 variability between PDX lines. Recall, this parameter was set to 100% for the flank 154 experiments. Intracranially, the model estimated that, on average, 20% of drug was 155 successfully delivered for GBM6, 48% for GBM12, and 100% for GBM39 (Fig 3B).

156
While GBM6 appeared to be the most sensitive of the PDX lines, the poor drug delivery 157 limited the response intracranially. Similarly, although GBM39 was estimated to be the 158 least sensitive of the PDX lines, the estimated high drug delivery resulted in a favorable 159 response to intracranial therapy. This suggests that drug delivery plays a major role in 160 treatment response. presented, by demonstrating the utility of the model through applying the parameter 171 estimation to serial BLI data from three PDX lines [7]. Broadly, we found the model 172 was able to fit data from multiple replicates of multiple cell lines. Through the key 173 parameter fittings of µ L and γ, the model provides insights into the different underlying 174 mechanisms driving the tumor growth behavior after therapy in the intracranial setting. 175 For GBM6, the model attributed the intracranial treatment failure to poor drug delivery, 176 presumably due to an intact BBB. While GBM12 was also fairly sensitive to the drug, 177 the model identified poor delivery as the major limitation, again presumably as a result 178 of the BBB. Unlike the other PDX lines, GBM39 responded well to treatment. The By using only serial BLI data, the model was able to make similar inferences about 181 the three PDX cell lines as the more rigorous results found in Marin et al. [6]. For 182 example, the model identified drug delivery as a key contributor to the outcomes of 183 GBM6 and GBM39. In fixed brain sections, elevated levels of fibrinogen are indicative 184 of BBB disruption and while GBM6 demonstrated no detectable fibrinogen 185 accumulation, GBM39 tumors had significant fibrinogen buildup near the region of the 186 tumor. Further, GBM6 tumors demonstrated minimal accumulation of drug while 187 GBM39 tumors exhibited drug accumulation near the region of the tumor [6].

188
The agreement, though not perfect, between these two approaches is very 189 encouraging. The BLI experiments used in this paper were conducted in parallel to the 190 experiments for Marin et al. [6]. While it was decided to publish Marin et al. first, the 191 results presented here were identified first [6]. We believe this demonstrates the fashion, streamlining initial cell line selection for more rigorous experiments [6]. In this 207 way, the model could reduce overall experimental costs and hopefully also decrease time 208 to meaningful answers for patients.