A Dynamical Systems Model for Combinatorial Cancer Therapy Enhances Oncolytic Adenovirus Efficacy by MEK-Inhibition

Oncolytic adenoviruses, such as ONYX-015, have been tested in clinical trials for currently untreatable tumors, but have yet to demonstrate adequate therapeutic efficacy. The extent to which viruses infect targeted cells determines the efficacy of this approach but many tumors down-regulate the Coxsackievirus and Adenovirus Receptor (CAR), rendering them less susceptible to infection. Disrupting MAPK pathway signaling by pharmacological inhibition of MEK up-regulates CAR expression, offering possible enhanced adenovirus infection. MEK inhibition, however, interferes with adenovirus replication due to resulting G1-phase cell cycle arrest. Therefore, enhanced efficacy will depend on treatment protocols that productively balance these competing effects. Predictive understanding of how to attain and enhance therapeutic efficacy of combinatorial treatment is difficult since the effects of MEK inhibitors, in conjunction with adenovirus/cell interactions, are complex nonlinear dynamic processes. We investigated combinatorial treatment strategies using a mathematical model that predicts the impact of MEK inhibition on tumor cell proliferation, ONYX-015 infection, and oncolysis. Specifically, we fit a nonlinear differential equation system to dedicated experimental data and analyzed the resulting simulations for favorable treatment strategies. Simulations predicted enhanced combinatorial therapy when both treatments were applied simultaneously; we successfully validated these predictions in an ensuing explicit test study. Further analysis revealed that a CAR-independent mechanism may be responsible for amplified virus production and cell death. We conclude that integrated computational and experimental analysis of combinatorial therapy provides a useful means to identify treatment/infection protocols that yield clinically significant oncolysis. Enhanced oncolytic therapy has the potential to dramatically improve non-surgical cancer treatment, especially in locally advanced or metastatic cases where treatment options remain limited.

We fixed the estimated values (Table S1) and fit MEK--inhibition induced G1 cell cycle phase arrest and release parameters (in red) to CI1040 data ( Figure 1b). Resulting simulations were within 6% of the corresponding data: dC/dt = σ·•C·•(1--P/sat) -arrest(t;t d ,t w )·•C + release·•C G1 dC G1 /dt = arrest(t;t d ,t w )·•C -release·•C G1 P = C + C G1. The parameter directing the rate at which cells release from G1 cell cycle arrest, release, was evaluated as both a constant parameter and a dynamic function. We retained the constant parameter form of release since resulting simulations offered greater precision (i.e., a lower SSE). Upon fitting σ, sat, arrest, and release, we quantified virus uptake and cell viability using replication--competent ONYX--015 to characterize the infection and cell death dynamics without use of MEK--inhibition ( Figure S2). As a result, we were able to use another reduced form of the model to fit β, δ, and corresponding delay terms (in red) to DMSO infection data: dC/dt = σ·•C·•(1--P/sat) -β(t;t i ,MOI)·•C dIC/dt = β(t;t i ,MOI)·•C -δ(t;t i ,MOI)·•IC P = C + IC We fit parameters to different versions of the model maintaining the form that included a unique delay term associated with infection, and a unique delay term associated with cell death. Parameter fitting yielded 4 unique values (per MOI) with corresponding kinetics within 3, 5, 6, 7, and 9--percent of the experimental data for MOIs of 0.1, 1, 2, 5, and 10, respectively ( Figure S2; Table S1). We fixed DMSO specific infection/viability parameters and fit the remaining virus uptake and cell viability parameters to experimental data that quantified the cellular response to replication competent ONYX--015 infection after 2 days of MEK--inhibitor pre--treatment ( Figure 1e). This data supported fitting the β T , β T·•G1 , δ T , and corresponding delay terms (in red) associated with the CI1040 portion of the model: T We fit delay parameters to different versions of the model: i. Unique delay term for δ T . ii.
Unique delay term for δ T , and a common delay term for β T and β T·•G1 . iv.
Unique delay terms for δ T , and β T . The model reflecting greatest precision with the data incorporated delays in infection and lysis when the population of cells was not in G1--arrest (v). In other words, model development suggested that cells in CI1040-mediated G1 arrest immediately infect and transition to the treated--infected state, IC T . Parameter fitting for this case yielded 5 unique values (per MOI) that were accurate to within 5, 8, 8, 7, and 6--percent of the experimental data for MOIs of 0.1, 1, 2, 5, and 10, respectively ( Figure S3). The resulting parameters (Table S1) were specific to cells that had been pre--treated with CI1040 for 2 days, followed by media change and infection.

Parameter Convergence
We assessed convergence of parameter estimation by evaluating the (W)SSE and resultant parameter values of 50 independent iterations of the optimization algorithm. If the algorithm returns consistent parameter estimates and consistent fitness (as denoted by the cost function), parameter identifiability issues are likely negligible since the optimization method is able to account for and estimate individual parameter values. Alternatively, if the algorithm returns consistent measures of fitness and variable parameter values, the model and algorithm may suffer from an inability to identify parameters independently. In the latter case, the cost function is relatively insensitive to corresponding parametric perturbations. A common cause of identifiability problems entails dependency among parameter values. A simple test involves assessing the correlation between resulting parameter values when such variance exists. Among the proliferation parameters, arrest and release exhibited clear mutual dependence with a correlation coefficient of 0.993 ( Figure S4). Their correlation was expected, however, since both parameters regulate the reversible G1 cell cycle arrest dynamics as a result of MEK--inhibitor. Similarly, sigma and sat exhibited a correlation coefficient of 0.821. Additional correlations were observed among certain infection and cell death kinetics in DMSO conditions. At higher MOI (i.e. 2, 5, and 10), the rates of infection and the rates of cell death are weakly correlated with coefficients between 0.201 and 0.353. The rates of infection or cell death are comparably correlated with their corresponding delay terms. Parameters that exhibited the strongest association involved the delay of infection and the rate of cell death (correlation coefficients between 0.604 and 0.796 across MOI). These observations are justified by the fact that cell death is an irreversible reaction that must take place once a cell is infected. At low MOI (i.e. 0.1 and 1), parametric correlations were inconsistent, which likely relates to the lack of experimental dynamics since few cells infect and lyse in low virus titer. CI1040 conditions also yielded inconsistent parameter correlations with respect to MOI. The one exception involved the rate of infection from the arrested state, β T·•G1 , and the delay of infection from the treated proliferating state (correlation coefficients ranged between 0.493 and 0.802). The variability of parametric correlation in the CI1040 case is likely related to the greater complexity of the model structure. Here, cells may infect into the treated state, IC T , via the proliferative state, C, or the arrested state, C G1 . Figure S4: Identifying correlation between parameters arrest and release. Fifty iterations of independent parameter estimation demonstrated a correlation between parameters arrest and release with a correlation coefficient of 0.993. Such a strong association was expected since arrest and release govern the only reversible reaction in the system of equations: G1 cell cycle arrest. Subplots on the left depict the histogram of the 50 resulting parameter values while the plot on the right depicts their correlation. The observed parametric correlations did not present adequate concern to motivate the derivation of a non-dimensional form of our nonlinear system of ODEs. Deliberate sequential parameter estimation required fitting a maximum of 5 parameters at any given time, limiting the uncertainty (or insensitivity) of the fitness function with respect to parameter values. To assess convergence, we computed each optimization algorithm 50 independent times. The resultant parameters exhibited negligible variation. Exceptions included both delay terms in the DMSO cases at low MOI, and parameters β T•G1 and delay δT in the CI1040 cases at low MOI. Given the limited experimental dynamics at MOIs of 0.1 and 1, we did not believe that reducing the number of fitted parameters would ensure greater precision or significantly affect the aims of our study.

Model Fitness
In using a maximum likelihood approach to estimate model parameters, we assumed that the residuals of our fitness function were normally distributed and zero mean centered. Therefore, model fitness is determined in part by whether this assumption holds. Figure S5 depicts the histogram of residuals resulting from each of the 50 independent iterations of parameter estimation. In the leftmost plot, residuals associated with sigma, sat, arrest, and release reflect a normal distribution at the 5% significance level with a mean of --0.90 (standard deviation of 6.94). In the middle plot, residuals associated with β, δ, and corresponding delay terms (DMSO treatment) reflect a normal distribution at the 5% significance level with a mean of 4.45 (standard deviation of 8.40). In the rightmost plot, residuals associated with β T , β T·•G1 , δ T , and corresponding delay terms (CI1040 treatment) reflect a normal distribution at the 5% significance level with a mean of 3.67 (standard deviation of 9.59). If all residuals are consolidated, they are normally distributed at the 5% significance level with a mean of 2.77 (standard deviation of 8.80). Residuals, x, meet the criteria for normal distribution as a function of MATLAB's t--test provided a mean, m: 'ttest(x,m) performs a t--test of the null hypothesis that data in the vector x are a random sample from a normal distribution with mean m and unknown variance'. It is important to note that residuals are not zero--mean centered and fail the test if the mean is not specified. This deviation may be associated with initial conditions since the time interval that governs the ODE solver may not be consistent for each experimental measurement. In particular, the first time point (or residual) may bias the SSE metric since it is anchored to the first experimental measure and subject to experimental variation and initial conditions unaccounted for in the system of equations. Provided a non--zero mean of residuals, we cannot ensure that parameter estimates are optimal. However, given that normalized mean values (with respect to standard deviations) are minor (--0.13, 0.53, and 0.38 respectively) and that mean--centered residuals are normally distributed, a maximum likelihood approach to estimate model parameters is suitable.

Interpolation Methods
By fitting a wide range of infection conditions we were able to establish a framework that allowed for interpolation of intermediate parameter values to account for a continuum of experimental and therapeutic possibilities. To assess the efficacy of MATLAB's 1D interpolation algorithms, we determined the total cost associated with interpolating known data points. More specifically, we omitted parameters corresponding to an MOI of 1, 2, or 5, interpolated the absent parameter value given the spectrum of known parameters, and calculated the SSE between the interpolated parameter simulation and the training data. The total SSE for these MOIs was determined for each relevant MATLAB interpolation method. Results demonstrated that cubic interpolation was most accurate for this system ( Figure S6).

Model Simulations
We interpolated resulting parameter values (Table S1) and used the model to predict the extent of cell death as a function of the time of CI1040 treatment initiation, time of ONYX--015 infection, and MOI. We employed an exhaustive search algorithm to simulate the effect of various treatment and infection protocols. This algorithm systematically evaluated every possible sequence combination of drug treatment and infection conditions (within a defined interval), with the exception of media change, t w , which was set to occur 2 days after treatment. We varied CI1040 treatment initiation between days 0--3 and infection between days 0--7. The multiplicity of infection was also varied between 0.1 and 10 ( Figure S7). We evaluated percent cell death on day 8 irrespective of the sequence protocol.

Minimum Information about a Flow Cytometry Experiment (MIFlowCyt) Outline
In this section, we outline how our flow cytometry experiments satisfy MIFlowCyt guidelines, as established by Lee et al. 1 . All FCS data files can be obtained by contacting Dr. Michael Korn. Figure 1A: