Leveraging Hypoxia-Activated Prodrugs to Prevent Drug Resistance in Solid Tumors

Experimental studies have shown that one key factor in driving the emergence of drug resistance in solid tumors is tumor hypoxia, which leads to the formation of localized environmental niches where drug-resistant cell populations can evolve and survive. Hypoxia-activated prodrugs (HAPs) are compounds designed to penetrate to hypoxic regions of a tumor and release cytotoxic or cytostatic agents; several of these HAPs are currently in clinical trial. However, preliminary results have not shown a survival benefit in several of these trials. We hypothesize that the efficacy of treatments involving these prodrugs depends heavily on identifying the correct treatment schedule, and that mathematical modeling can be used to help design potential therapeutic strategies combining HAPs with standard therapies to achieve long-term tumor control or eradication. We develop this framework in the specific context of EGFR-driven non-small cell lung cancer, which is commonly treated with the tyrosine kinase inhibitor erlotinib. We develop a stochastic mathematical model, parametrized using clinical and experimental data, to explore a spectrum of treatment regimens combining a HAP, evofosfamide, with erlotinib. We design combination toxicity constraint models and optimize treatment strategies over the space of tolerated schedules to identify specific combination schedules that lead to optimal tumor control. We find that (i) combining these therapies delays resistance longer than any monotherapy schedule with either evofosfamide or erlotinib alone, (ii) sequentially alternating single doses of each drug leads to minimal tumor burden and maximal reduction in probability of developing resistance, and (iii) strategies minimizing the length of time after an evofosfamide dose and before erlotinib confer further benefits in reduction of tumor burden. These results provide insights into how hypoxia-activated prodrugs may be used to enhance therapeutic effectiveness in the clinic.


Derivation of erlotinib plasma concentration function
We define a function ρ 1 (t) describing plasma concentration over time after one dose of D 1 mg erlotinib using multiple sets of pharmacokinetic data. Assume the concentration of erlotinib decays exponentially as a function of time and let ρ 1 (t) = c max 1 e −k1t describe this exponential decay, where c max 1 is the maximum plasma concentration achieved after one dose of erlotinib and k 1 is the rate of decay. Using available pharmacokinetic data, we acquire pairs of best-fit parameters for this function at a variety of doses between 100 and 300 mg. The values obtained for k 1 are averaged to produce k 1 = 0.0307, and a linear function c max 1 (D 1 ) = 7D 1 + 365 is defined to reflect the relationship between dose and maximum plasma 1 concentration. Hence we define erlotinib plasma concentration, given an initial dose D 1 , as a function of time: ρ 1 (t) = (7D 1 + 365)e −0.0307t .

Derivation of birth and death rates due to evofosfamide
Here we define functions for the birth and death rates of sensitive and resistant cells in terms of evofosfamide concentration. We start by defining the percent of viable cells V i (C 2 ) in compartment i as a function of evofosfamide concentration. Suppose that for a fixed concentration of oxygen, cell viability behaves as a logistic function with respect to drug concentration. Let V (C 2 ) = 100 − A 1+e −α/β + A 1+e (C 2 −α)/β , with A, α, and β unknown parameters, to reflect the fact that cells are 100% viable in the absence of drug (V (0) = 100). We fit this logistic function to multiple experimentally-derived data sets for cell viability, each at a different fixed oxygen concentration between 0% and 21%, to obtain values for A, α, and β. The estimates of A are averaged to yield A = 88.5983, and the relationship between α and β for compartment i and x i , the oxygen concentration in compartment i, is described as follows: β i = (6.389 · 10 −8 )x 3 i − (9.284 · 10 −7 )x 2 i + (3.87 · 10 −6 )x i + 3.792 · 10 −8 x i < 10% (6.366 · 10 −7 )x i − 3.421 · 10 −6 x i ≥ 10%.
Therefore, for compartment i, the percent of viable cells at an evofosfamide concentration of C 2 is expressed as where A, α i , and β i are all defined as above. Fig. 1 shows a plot of this function together with the corresponding data points at the seven fixed oxygen concentrations examined in the cell viability experiments.
We now use experimental cell viability data together with the cell viability function to derive expressions N X,i (C 2 ) and N Y,i (C 2 ) for net growth rates of sensitive and resistant cells in compartment i as functions of evofosfamide concentration C 2 . Let the fraction of viable cells in compartment i as a function of evofosfamide concentration be denoted by v i (C 2 ) = 1 100 V i (C 2 ). Since we assume evofosfamide has the same effect on the sensitive and resistant cells, we need only define one function for net growth rate in compartment i in terms of evofosfamide concentration which we can then apply to both the sensitive and resistant cell populations. Let N (x, d) represent the net growth rate of cells in an environment with oxygen concentration x (written as a percentage) and evofosfamide concentration d. We define this function by writing out explicit expressions for the number of cells in the treated and control groups at certain points in time according to the details of the cell viability experiments and comparing this with the cell viability function.
Let group A represent the control group in the cell viability experiments and let group B denote the group of cells treated with evofosfamide concentration C 2 at an We use the detailed process followed in the cell viability experiments to write explicit expressions for the number of cells in groups A and B at certain points in time. At time t = 0, there are n A (0) = n B (0) cells in each group. According to the cell viability data, exponentially growing cells were seeded 24 hours before the addition of test compounds. We express the number of cells present in each group at t = 24 by We proceed by relating the above equations to cell viability. Setting t = 98 in Eq. (1) yields By substituting the equations in (2) into Eq. (3) and simplifying, keeping in mind that n A (0) = n B (0), we obtain Taking the natural log of both sides, dividing by 2, and solving for N (x i , C 2 ) yields We can rewrite this, applying the above results to both sensitive and resistant cells, to express net growth rates in compartment i as functions of evofosfamide concentration C 2 : Note that in the above equations, N X,i (0) and N Y,i (0) are simply the net growth rates of sensitive and resistant cells in the absence of drug. In other words, where the superscript c denotes the control growth rates in the absence of drug, as defined in the first section of this document.
However, having explicitly defined functions for net growth rates of the sensitive and resistant cell populations is not quite enough; we must define their respective birth and death rates as well. Because evofosfamide is a cytotoxic drug, we assume it has no effect on birth rates, but rather increases death rates. Hence the birth rates of cells under any concentration of evofosfamide are equal to those in the absence of drug. Thus we define birth rates of sensitive and resistant cells in compartment i due to evofosfamide by Using the fact that net growth rates are differences between birth and death rates along with a series of substitutions of the above equations followed by algebraic simplification yields the following functions for death rates of sensitive and resistant cells in compartment i:

Derivation of evofosfamide plasma concentration function
Using pharmacokinetic data from a phase 1 clinical trial, we define a function ρ 2 (t) for plasma concentration over time after one dose of D 2 mg/m 2 evofosfamide. Unlike with erlotinib, we cannot simply approximate the plasma concentration of evofosfamide using only exponential decay because doing this would force the maximum concentration to be artificially high since evofosfamide decays so quickly with respect to time. This does not present a problem in the erlotinib plasma concentration function because the decay rate of erlotinib is much lower than that of evofosfamide. From available clinical data, we estimate the time at which maximum plasma concentration is achieved to be t = 1 2 (30 minutes after the initial dose is administered). We approximate the plasma concentration for t ≥ 1 2 as exponential decay using ρ 2 (t) = c max 2 e −k2(t−1/2) , where c max 2 is the maximum plasma concentration achieved after one dose and k 2 is the decay rate of evofosfamide. We acquire pairs of best-fit parameters for this function at a variety of doses between 7.5 and 940 mg/m 2 using the clinical trial data. Cubic functions are defined for both c max 2 and k 2 to describe the relationship between dose D 2 and maximum concentration and decay rate as follows: To estimate the plasma concentration for t < 1 2 , we define a linear function connecting (0, 0) with 1 2 , c max 2 (D 2 ) . We can now define evofosfamide plasma concentration after one dose as a function of time as follows: (D 2 ) and k 2 (D 2 ) are each defined as in the above equations.

Definition of one cycle of treatment in a combination dosing schedule
Here we describe in detail the process by which a single cycle of treatment in a combination dosing schedule is defined. We begin by calculating the length L of the cycle, given n. Since n is the number of evofosfamide doses given in 3 weeks and there are 504 hours in 3 weeks, L = 504/n. Next we calculate the maximum tolerated dose of evofosfamide, given n, from the evofosfamide toxicity constraint curve. This dose of evofosfamide is administered at time t = L − 24 for Classes 1 and 2, and t = L − 6 for Class 3. The remaining time in the cycle, starting at t = 0 and going up to the time of the evofosfamide infusion, is then filled with the fixed erlotinib dosing schedule associated with the given optimization class. If necessary, doses of erlotinib are then removed, beginning with the dose immediately before evofosfamide is administered and working backward, until there is enough time between the last dose of erlotinib and the first dose of evofosfamide for the erlotinib plasma concentration to fall below 2.357 µM so as to not violate the combination toxicity constraint.

Explanation of the presence of two local minima in the optimization results
One interesting characteristic of the means of the sensitive and resistant cells in the optimization results is the presence of two local minima. Here we show that this phenomenon is due to the variability in timing between the last dose of erlotinib and the first dose of evofosfamide in each cycle.
To eliminate the variation in timing between doses, we re-define all dosing schedules in Class 1 so that either erlotinib or evofosfamide is administered every day with no additional doses in between. After calculating the length L of one cycle, we round L − 24 to the nearest 24 hours and define that to be the time of the evofosfamide infusion. Then, ignoring combination toxicity constraints solely for the purpose of this argument, we schedule erlotinib doses every 24 hours prior to the evofosfamide infusion. After re-defining dosing schedules in this manner, there no longer exists variability in the time between drug administrations since there is always a 24-hour waiting period between any two doses.
Originally, because the n evofosfamide doses were evenly divided into every three weeks, the end of treatment (at nine weeks) was also the end of a treatment cycle for every dosing schedule. However, that is no longer the case for these new daily dosing schedules since doses are forced to occur every 24 hours and only every 24 hours. Thus there are now two options for what to consider as the end of treatment. We may define the end of treatment to be exactly nine weeks from the treatment start date, but this introduces variability in the position within the cycle at which the resistance dynamics are calculated. The other option is to define the end of treatment to be the end of the cycle closest to nine weeks, but this will be at different times for different dosing schedules.
We examine both of these possibilities by re-calculating the means of the sensitive and resistant cells at both of these times for Class 1 under these new daily dosing schedules, shown in Fig. 2A and Fig. 2B, respectively. For the sake of comparison, the optimization results from the original dosing schedules in Class 1 are shown in blue. For each daily dosing schedule, the means of the sensitive and resistant cells are calculated at two different times: nine weeks from the start of treatment (yellow) and the end of the cycle closest to nine weeks after treatment begins (red). During the creation of the daily dosing schedules, the rounding of L − 24 to the nearest 24 hours leads to multiple dosing schedules with the same dose timing but different doses. To simplify the analysis, we only include one dosing schedule with each particular schedule of doses in Fig. 2, chosen based on which n yields the value of L − 24 closest to the common rounded value.
Let us now compare the results from the original dosing schedules with those from the daily dosing schedules. The elimination of the variability in time between doses causes the first local minimum to become less pronounced for both the mean of the sensitive cells and the mean of the resistant cells, regardless of how the end of treatment is defined. Although the first local minimum still exists, its existence can now be easily explained in every case based on the different types of variation that have been introduced.
When the end of treatment is defined to be nine weeks after the treatment start date, this introduces variability in the relative position within the cycle at which the evolutionary dynamics are calculated. Specifically, the point at which this calculation occurs during treatment for n = 3 is at the very end of a cycle (after a large decline in the cancer cell population due to evofosfamide), whereas the point at which this calculation occurs for n = 2 and n = 4 is right before a dose of evofosfamide. This explains the existence of a local minimum at n = 3 for both means when the end of treatment is defined to be nine weeks after the start of treatment.
On the other hand, when the end of treatment is defined to be the end of a cycle, this introduces variability in the actual time at which the recurrence dynamics are calculated, and hence variability in the total treatment time. Specifically, the total treatment time for n = 2 is quite high relative to other values of n. Since the sensitive cell population is in a constant state of decline during treatment, this explains the existence of a local minimum at n = 2 for the mean of the sensitive cell population when the end of treatment is defined as the end of a cycle. Additionally, the total treatment time for n = 3 is quite low relative to other values of n. Since the resistant cell population is steadily growing by the end of treatment, this explains the existence of a local minimum at n = 3 for the mean of the resistant cells when the end of treatment is defined as the end of a cycle.
Note that for the mean of the sensitive cells, by following the curve for which the end of treatment is defined to be at nine weeks up to n = 3 and then the curve for which the end of treatment is defined to be the end of a cycle, the first local minimum is completely eliminated. This is a reasonable comparison to make because for n = 3, nine weeks after the start of treatment is the end of a treatment cycle. Together, all of the above evidence shows that the first local minimum is due to the variability in timing between the last dose of erlotinib and the first dose of evofosfamide in each cycle. treatment for a tumor initially consisting of 1.6 · 10 6 sensitive cells are plotted for the original dosing schedules with variation in time between doses from Class 1 (blue) as well as the newly-defined daily dosing schedules, which eliminate this variation. For each daily dosing schedule, these means are calculated at two different times: nine weeks from the start of treatment (yellow) and the end of the cycle closest to nine weeks after treatment begins (red).

Impact of migration
To investigate the impact of migration on our model predictions, we relaxed the assumption that the evolutionary dynamics within each microenvironmental compartment are independent by introducing a migration term that allows cells to move between compartments. In addition to proliferating and dying, suppose cells in compartment i migrate to compartment j with rate ρ(i, j). Assume that cells in compartment i can only migrate to compartment i + 1 or i − 1 at a given point in time. Then the means of the sensitive and resistant cells in each compartment are governed by a system of 64 ODE's. For i = 1, 2, ..., 32, we have In the above equations, let ρ(1, 0) = ρ(0, 1) = 0 and ρ(32, 33) = ρ(33, 32) = 0 to account for the fact that there are only 32 compartments. We defined the migration rate ρ(i, j) using three different methods. Let the rates defined by these methods be denoted as ρ A (i, j), ρ B (i, j), and ρ C (i, j). For each compartment, by using the decay rate of oxygen away from the blood vessel together with the range of oxygen concentrations in that compartment, we estimated the range of distances from the blood vessel that correspond to the given compartment. We then used the minimum and maximum distances for each compartment to estimate a measurement of "width" -call this w i for compartment i. To estimate w i for the compartment furthest from the nearest blood vessel, we used the fact that the mean inter-capillary distance in tumor tissue is approximately 300 µm, so the point furthest from the nearest blood vessel in this compartment is at a distance of approximately 150 µm [2,3]. We assume that cells travel at a rate of r = 0.3 µm/min [4]. The first method we used to define ρ A (i, j) assumes that cells only migrate toward the nearest blood vessel. So the migration rate is given by The second method we used assumes that cells only migrate away from the nearest blood vessel. Then the migration rate is given by For the third method, we allowed cells to migrate both toward and away from the nearest blood vessel with a constant migration rate of ρ C (i, j) = 0.1. This is a reasonable estimate since the order of magnitude of r is −1 and for most compartments, w i has order of magnitude 0.
To examine the impact of migration on the mean tumor size, we calculated mean tumor size over time for two different combination dosing schedules using migration defined in each of the three methods described above as well as without migration. Fig. 3A shows the results of this calculation for a dosing schedule that combines a standard 150 mg daily erlotinib schedule with a standard 575 mg/m 2 weekly evofosfamide schedule. Fig. 3B shows the results for the combination dosing schedule given by n = 14 in Class 3, which was found to be the optimal dosing schedule to minimize mean tumor size in the previous section.
We find that migration does have an impact on our model predictions with regards to mean tumor size. Since the majority of the cancer cells are initially in hypoxic regions of the tumor microenvironment, allowing migration between compartments governed by a constant migration rate leads to a more well-mixed population. Therefore, with a constant migration rate incorporated into the model, erlotinib has more of an effect on the population while evofosfamide has less of an effect since the cells are more evenly distributed throughout the tumor microenvironment. This can be seen in Fig. 3A and Fig. 3B. This same phenomenon can be observed in the case where cells only migrate toward the blood vessel, although it is more pronounced in this situation since cells are only allowed to move away from the hypoxic regions of the tumor microenvironment. On the other hand, by incorporating migration such that cells may only move away from the blood vessel, the observed shift in mean tumor size is likely due to a decreased response to erlotinib and an increased response to evofosfamide due to the fact that the cancer cells are even more concentrated in hypoxic regions.
These results indicate that unidirectional migration of cell populations either towards or away from blood vessels have more of an impact on the model-predicted mean tumor evolution than bidirectional migration rates. However, the migration rates chosen here fell on the higher end of literature-reported values, and it is likely that migration patterns found in vivo are multidirectional, heterogeneous and dependent upon cellular density in local and surrounding compartments. In addition, it is possible that sensitive and resistant cells may have different migration speeds or capabilities. To better quantify the impact of migration on the evolutionary dynamics of the tumor population, further experimental investigation is needed to quantify migration rates and develop a more refined mathematical model of local migration between compartments.

Supplementary matrices
Here we show the matrices SB, SD, RB, and RD used to calculate the birth and death rate functions of erlotinib-sensitive and erlotinib-resistant cells, defined in the first section of this document.