Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution

The success of targeted cancer therapy is limited by drug resistance that can result from tumor genetic heterogeneity. The current approach to address resistance typically involves initiating a new treatment after clinical/radiographic disease progression, ultimately resulting in futility in most patients. Towards a potential alternative solution, we developed a novel computational framework that uses human cancer profiling data to systematically identify dynamic, pre-emptive, and sometimes non-intuitive treatment strategies that can better control tumors in real-time. By studying lung adenocarcinoma clinical specimens and preclinical models, our computational analyses revealed that the best anti-cancer strategies addressed existing resistant subpopulations as they emerged dynamically during treatment. In some cases, the best computed treatment strategy used unconventional therapy switching while the bulk tumor was responding, a prediction we confirmed in vitro. The new framework presented here could guide the principled implementation of dynamic molecular monitoring and treatment strategies to improve cancer control.

1 Mathematical Methods

Evolutionary Dynamics Model of NSCLC
The quasispecies model [1] was originally developed to describe the dynamics of populations of self replicating macromolecules undergoing mutation and selection. We choose this model for its relative simplicity and its ability to capture the salient features of the evolutionary dynamics of a simplified generic disease model. The following adaptation incorporates the effects of small molecule inhibitors and describes the growth, mutation and evolution of non small cell lung adenocarcinoma populations: where x i ∈ R + is the concentration of a NSCLC subpopulation i, k ∈ R + is a small molecule inhibitor concentration (assumed to remain at constant concentrations throughout), r i is the growth rate for each cell x i , and q ik is the probability that cell k mutates to cell i (note that q ii is the probability of no mutation occurring). Finally, the function Ψ i ( k ) represents the pharmacodynamics of individual drugs k or of individual EGFR TKIs (erlotinib or afatinib) in combination with fixed concentrations other small molecule inhibitors used in this study (0.5 µM crizotinib, 0.5 µM trametinib or 5 µM vemurafenib) with respect to the i-th NSCLC cell type, namely: where k ∈ R + is the drug concentration, γ ik ∈ R + is the saturation coefficient, K ik ∈ R + is the dissociation constant, n k ∈ R + is the Hill coefficient. Equation S2 has previously been described in [2,3]. When k = 0, ∀k ∈ {1, ..., m}, the dynamics are unstable.

A control theoretic algorithm for designing treatment strategies
To design treatment strategies that best minimize tumor size and control its evolution over time, we combine both a greedy algorithm and receding horizon control approach. We introduce some notation, cost function definitions and specify our algorithm.

Cost functions
To measure the effectiveness of a given treatment strategy over time, we define the average cost function. For a given treatment strategy k applied to Equation (S1), we rewrite the dynamics of the entire system (i.e., for all cells) asẋ where A ∈ R n×n is a matrix that represents the growth, mutation and drug dynamics for treatment strategy k , for n cell subpopulations.
The average cost C r for a time horizon N , allowable switching period τ and time intervals of the form [kτ, (k + 1)τ ] for k = {0, .., N/τ − 1} is given by where 1 T is the n × 1-dimensional vector of ones and x(t) is the solution to Equation (S3). Equation (S4) simplifies to The final cost C f for an inital tumor population x(0) and a sequence of drugs { (k)} that define a switching therapy over a time horizon N is defined as

Algorithm
Our algorithm is defined as follows. Given an initial tumor population, denoted by x 0 , a time horizon N and an allowable switching period τ , we perform the following computations to determine a candidate treatment strategy: Algorithm 1 Treatment strategy synthesis 1. Initialization: Set k = 0 and x(0) = x 0 .
that define a switching therapy.
The resulting switching therapy { (k)} is then applied until the next biopsy can be taken, giving a new tumor cell population measurement, at which point the algorithm is repeated. In particular, it is important that the horizon N be chosen to be longer than expected periods between biopsies.

Derivation of dynamical system parameters
Growth and Mutation Rates. We model the growth of NSCLC cell population x i by the following ordinary differential equation (ODE):ẋ where r i is the growth rate per day, andẋ i denotes the derivative with respect to time of the tumor cell population x i . Note that we assume that no mutations occur over the time-frame considered, allowing us to set q ii = 1 and q ij = 0 in the dynamic model (S1), resulting in (S7). Given an initial population x i (0), the population x i (t) on day t can be obtained by solving ODE (S7), and is specified by the following expression Given a set of N experimental data points e i (0), e i (t 1 ), . . . , e i (t N ), we fit these points to an exponential function of the form (S8), with x i (0) = e i (0) to obtain an experimentally derived value for the growth rate r i of tumor cell population x i .
We take the DNA mutation rate to be 1e −9 mutation/base pair/cell division. We assume that mutations occur unidirectionally from EGFR L858R parental cells to EGFR L858R,T790M , EGFR L858R , BRAF V600E or EGFR L858R,T790M BRAF V600E , HGF-/+. For a NSCLC cell population with growth rate r i , the corresponding doubling time td i (cell division per day) is The mutation rate in units of mutation/base pair/day for an NSCLC cell population with doubling time td i is 1e −9 · td i −1 . The rate of mutation to one particular base pair/day is then approximated by 1e −9 · g −1 s · td i −1 , where g s = 3e 9 is the size of the human genome in base pairs.
Drug Effect Rates and Hill Functions. We model the change in a tumor cell population x i under a treatment j of concentration with the following ordinary differential equation (ODE): where r i is the growth rate per day derived in the previous section and f j i ( ) is a function mapping the treatment j at concentration to a drug effect rate per day. We again assume that no mutation occurs over the time-frame considered, allowing us to set the mutation rates q ii = 1 and q ij = 0 in the model (S1), resulting in (S10).
Similar to the previous section, given an initial population x i (0), the population x i (t) on day t can be obtained by solving ODE (S10), and is specified by the following expression (S11) We model the map f j i ( ) as a modified function of the form where γ j,i n j,i and K j,i are the saturation parameter, Hill function coefficient and binding reaction dissociation constant for drug j applied to cell x i . Our goal is to obtain values for these three parameters using experimental data measuring cell viability under varying concentrations of drug j. In particular, given experimentally obtained data pairs of the form , y i,j, (1), where y i,j, (1) is the ratio of the tumor cell population x i treated with concentration of drug j at day 1 to the tumor cell population x i treated with no drug at day 1. Letting x i denote the treated tumor population and x ctrl denote the untreated control tumor population, it follows that y i,j, can be written as where the first equality follows from the definition of y i,j, (1), the second from applying equations (S11) and (S8) to x i (1) and x ctrl i (1) respectively, and the third from canceling like terms. It follows that the experimentally derived values of f j i ( ) are given by Solving this equation for each experimentally tested concentration , we obtain a set of points { , f j i ( )} that can be used to derive the parameters γ j,i n j,i and K j,i via curve fitting. In order to avoid overfitting, we set γ j,i = max f j i ( ), i.e., we force the modified Hill function to saturate at the maximal experimentally observed rate. Although this approach can be conservative in modeling the drug effect rate of high concentrations of drugs, we note that the the maximal dose tested is chosen to be significantly higher than the maximum tolerated doses, and hence we do not expect this saturation to affect the accuracy of our model at clinically relevant doses.

Evolutionary stability measured by maximum eigenvalues
where x ∈ R n is a vector of concentrations of n NSCLC subpopulations,ẋ ∈ R n is their rate of change over time, A ∈ R n×n is a matrix that represents the growth and mutation dynamics and D ∈ R n×n is a diagonal matrix that represents the corresponding drug effect of one constant drug treatment on the rate of change of NSCLC cells. If all eigenvalues are negative then Equation (S15) is said to be stable. In the case of NSCLC evolutionary dynamics corresponding to Equation (1), stability refers to tumor reduction, and instability refers to tumor progression. In section 3.1, we made the assumption that mutation rates are one directional, hence the A matrix in Equation (S15) is lower triangular and the eigenvalues of A − D are exactly equal to its diagonal entries. For each NSCLC subpopulation, we take the maximum eigenvalue for each evolutionary branch downstream of the population and define this as evolutionary stability. This maximum eigenvalue represents the worst case stability if the particular population is present upon treatment initiation -a positive maximum eigenvalue indicates that the presence of the cell subpopulation in the tumor upon initiation of treatment is likely to cause therapeutic failure. A negative maximum eigenvalue indicates that the presence of the particular subpopulation will not outgrow or evolve in the presence of therapy.

Robustness analysis
Sensitivity to drug perturbations. To analyze the effect of dose reductions on the robustness of constant and switching treatment strategies, we perturbed the drug concentrations and calculated the ratio of final cost and initial cost (Figures (S9)) . We rewrite Equation (S1) for one cell x i and one drug j to illustrate how a drug perturbation δ ∈ R [0,1] is modeled: The fold change F C f in total population from day 0 to day N for a sequence of drugs defining a switching strategy over a time horizon N , and initial tumor population is effective for NSCLC populations for the duration of the time horizon N , F C f > 1 indicates progression.

Implementation
The evolutionary dynamics model and simulations were implemented using python, scipy and numpy (versions 3.5.1, 0.17.0, 1.9.3) and pandas version 0.17.0 was used for data parsing. Data fitting for experimentally derived cell growth and drug dose response data was performed with Matlab version 8.3.0.532 using the non linear least squares method.    Figure S4: Western blot analysis of cell lysates obtained from 11-18 cell line, treated with drugs and/or HGF as indicated, and probed for the indicated proteins. Corresponding full length gels are shown in S13, S14 and S15.  Figure S5: A) Experimentally derived dose response curves for afatinib in combination with 5 µM vemurafenib, 0.5 µM trametinib and 0.5 µM crizotinib for 11-18 EGFR L858R , 11-18 EGFR L858R BRAF V600E , H1975 EGFR L858R,T790M H1975 EGFR L858R,T790M BRAF V600E cell lines, and either 0 or 50 ng/ml human growth factor (HGF) and fit with γ [ ] n [ ] n +K n where γ is the maximum inhibition, [ ] is the EGFR TKI concentration, n is the Hill coefficient and K is the half maximal inhibitory concentration (IC50). (B) Western blot analysis of cell lysates obtained from H1975 cell lines, treated with drugs and/or HGF as indicated, and probed for the indicated proteins. Corresponding full length gel is shown in S16.   Figure S7: Simulations of the NSCLC model for the optimal 30 day constant combinations found by Algorithm (4) with 0.5 µM afatinib or 1.5 µM erlotinib with either 0.5 µM trametinib, 0.5 µM crizotinib or 5 µM vemurafenib for the relatively low (A) initial tumor heterogeneity or with (B) high initial tumor heterogeneity.        Figure S4C (note, blots were cut prior to incubation with primary antibody).  Figure S16: Uncropped western blots for Figure S5B (note, blots were cut prior to incubation with primary antibody).       Table 7: Differential equation parameters as derived using Equation (S12), corresponding to experimentally derived dose response curves of erlotinib in combination with either 0.5 µM crizotinib, 0.5 µM trametinib or 5 µM vemurafenib for parental and engineered 11-18 EGFR L858R -positive lung adenocarcinoma cells.  Table 8: Differential equation parameters derived using Equation (S12), corresponding to experimentally derived dose response curves of afatinib in combination with either 0.5 µM crizotinib, 0.5 µM trametinib or 5 µM vemurafenib for parental and engineered 11-18 EGFR L858R -positive lung adenocarcinoma cells.