Computational and Mathematical Methods in Medicine

Volume 2018, Article ID 8760371, 14 pages

https://doi.org/10.1155/2018/8760371

## Developing a Minimally Structured Mathematical Model of Cancer Treatment with Oncolytic Viruses and Dendritic Cell Injections

^{1}Department of Mathematics & Statistics, The College of New Jersey, Ewing, New Jersey, USA^{2}Department of Mathematics and Computer Science, University of Richmond, Richmond, Virginia, USA

Correspondence should be addressed to Jana L. Gevertz; ude.jnct@ztreveg

Received 22 May 2018; Accepted 6 September 2018; Published 30 October 2018

Academic Editor: Filippo Castiglione

Copyright © 2018 Jana L. Gevertz and Joanna R. Wares. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Mathematical models of biological systems must strike a balance between being sufficiently complex to capture important biological features, while being simple enough that they remain tractable through analysis or simulation. In this work, we rigorously explore how to balance these competing interests when modeling murine melanoma treatment with oncolytic viruses and dendritic cell injections. Previously, we developed a system of six ordinary differential equations containing fourteen parameters that well describes experimental data on the efficacy of these treatments. Here, we explore whether this previously developed model is the minimal model needed to accurately describe the data. Using a variety of techniques, including sensitivity analyses and a parameter sloppiness analysis, we find that our model can be reduced by one variable and three parameters and still give excellent fits to the data. We also argue that our model is not too simple to capture the dynamics of the data, and that the original and minimal models make similar predictions about the efficacy and robustness of protocols not considered in experiments. Reducing the model to its minimal form allows us to increase the tractability of the system in the face of parametric uncertainty.

#### 1. Introduction

For many solid tumors, the most utilized cancer treatments are surgery, chemotherapy, and radiotherapy [1]. While this approach can effectively reduce tumor burden in the short term, long-term recurrence is the norm. This failure of conventional treatment modalities has spurred efforts to design novel cancer therapeutics.

One emerging treatment modality is oncolytic virotherapy. This technique involves targeting cancer cells using oncolytic viruses (OVs), standard viruses genetically engineered to replicate selectively in cancer cells. OV replication within a cancer cell creates a viral burden too large for the cell to support, which eventually causes the infected cancer cell to lyse [2]. The OVs that get released from the lysed cancer cell are then free to infect additional cancer cells. The lysing effects of OVs, while powerful, are also transient. In a clinical setting, OVs have generally proven to be insufficient to fully and permanently eradicate a solid tumor mass [2].

Another treatment modality receiving attention is gene therapy, defined as the introduction of genes of interest into cancer cells with therapeutic intent [3]. Gene therapy has been attempted using genes that mediate the release of cytokines, tumor suppressor genes, and apoptosis-related genes, to name a few. Independent of the gene of interest, this modality requires a vector that can efficiently deliver, and uniformly distribute, the gene product to solid tumors [4]. Oncolytic viruses can be used as such a gene-delivery vector.

The ability of oncolytic virotherapy to induce tumor cell lysis and to stimulate an antitumor immune response in a preclinical setting has led to a number of clinical trials for different tumor types. As of 2016, twenty different viruses have been studied as candidates for oncolytic virotherapy, and new candidate viruses continue to be studied [5]. In 2015, the U.S. Food and Drug Administration approved the first oncolytic virus therapy, T-VEC (Imlygic™), for the treatment of advanced melanoma [5]. A number of clinical trials using this and other oncolytic viruses are underway in both solid and liquid tumors, and while they tend to be incredibly well tolerated, to date efficacy has been inferior to other available therapies [5, 6]. That said, researchers continue to actively study the anticancer effects of oncolytic virotherapy, both as a drug to be used in combination with other modalities, and as a potential cancer vaccine [5, 6].

In the present work, we focus on preclinical data and model the use of OVs (the adenovirus (Ad) in particular) to deliver genes that boost the immune system’s ability to identify, target, and kill cancer cells. The transgenes of interest in this study are 4-1BB ligand (4-1BBL) and interleukin (IL)-12. 4-1BB is a costimulatory member of the tumor-necrosis factor receptor superfamily that is expressed on activated CD4+ and CD8+ T cells [7]. The binding of 4-1BB to its ligand, 4-1BBL, promotes the outgrowth of type-1 T helper cells and cytolytic effector T cells [4]. IL-12 is a cytokine that strongly stimulates the differentiation of naïve CD4+ T cells to type-1 T helper cells. IL-12 has been determined to be one of the most efficient antitumor cytokines in experimental animal models [8]. For conciseness, we refer to oncolytic adenoviruses concurrently acting as a vector for both the 4-1BBL and IL-12 transgenes as Ad/4-1BBL/IL-12.

Recent preclinical work of Huang et al. has shown that Ad/4-1BBL/IL-12 can cause tumor regression in a mouse model of melanoma [4]. This debulking is a consequence of both tumor cell lysis, as well as immune system stimulation resulting from the local release of 4-1BBL and IL-12. Ad/4-1BBL/IL-12 can also be combined with intratumorally injected dendritic cell (DC) vaccines, resulting in a greater antitumor response than elicited by either treatment alone [4]. DCs are immune cells that present antigens to other cells of the immune system. Antigen presentation triggers an adaptive immune response that results in the immune system actively seeking out cells expressing the presented cancer antigen [9]. Huang et al. developed DC vaccines by harvesting DCs from the bone marrow of tumor-bearing mice, and exposing them ex vivo to tumor-associated antigens until maturation [4].

Given the combined effectiveness of Ad/4-1BBL/IL-12 and DC injections, it is natural to ask in what order, and at what dose, one should administer these therapeutics to elicit the maximal antitumor response. This is an experimentally time-consuming and costly question to address. Mathematical-modeling techniques can help answer questions about complex biological systems without the associated experimental costs [10]. Differential equation models (frequently paired with experimental data) have been successfully used to improve treatment protocols involving oncolytic viruses [11–16]. Previously, we hierarchically developed and fit a mathematical model to the experimental data in Huang et al. [4]. This system of ordinary differential equations, involving six variables and fourteen non-initial-condition parameters, was shown to well describe the dynamics of OVs enhanced with one or more immunostimulatory molecules (4-1BBL, IL-12, or both), DC injections, and DC injections coupled with Ad/4-1BBL/IL-12 [17–19]. Note that for each treatment protocol, the model was fit to the *average* tumor volume data (averaged over the 8-9 mice in the treatment cohort, with some consideration of the standard deviation in the data).

Using the best-fit parameters obtained from the hierarchical fitting to the average, we previously discovered that administering three doses of OVs followed by three doses of DCs (OV-OV-OV-DC-DC-DC) is the optimal drug ordering [18] when constrained to considering the drug dosing and spacing used by Huang et al. [4]. Further analysis, however, led us to doubt the robustness of this prediction. First, we found that the doses of Ad/4-1BBL/IL-12 and DC used in Huang et al.’s experiments [4] were near a bifurcation point; that is, slightly altering the dose or sequence could drastically change the efficacy of the protocol [18]. A full-scale robustness analysis of optimal protocols using the Virtual Expansion of Populations for Analyzing Robustness of Therapies (VEPART) procedure confirmed that our originally predicted optimal strategy is fragile [19].

The fragile nature of the optimal protocol raises doubts about whether it will actually be effective in individual mice in the experimental population, as individual mice generally have dynamics that deviate from the average. And, in the era of personalized medicine [20, 21], there is a strong emphasis on tailoring treatment protocols to individual patients. However, the difficulty of collecting and analyzing patient-specific data, especially in the face of intratumor and temporal heterogeneity, makes personalizing therapy a real challenge [22]. Even individualizing a mathematical model is a highly nontrivial task, as it often requires finding the best-fit parameters in a high-dimensional parameter space, given a very limited amount of data about the individual. Therefore, before we can explore the challenging question of personalizing therapy using our model, we ask the following questions: does our mathematical model require all six variables and fourteen parameters to adequately describe the data? Or, can we simplify the structure of the model (number of variables and parameters) in order to make the model more amenable to personalization while retaining the goodness of fit to the experimental data? These are the questions that we answer in this work.

This paper is organized as follows. First, we briefly introduce the reader to the previously developed mathematical model of tumor growth and treatment with Ad/4-1BBL/IL-12 and DCs [17, 18]. Second, we introduce a collection of methods that we employ to test whether our original model is minimal in its structure. Third, we argue that the original model is not of minimal structure, supported by evidence from parameter 95% credible intervals, local and global sensitivity analyses of parameters, and a soft/stiff parameter analysis. This leads us to propose a minimal model that contains five variables (one less than the original model) and eleven non-initial condition parameters (three less than the original model). Next, we show that this minimal model fits the experimental data as well as the original model, and we further argue that the model is not too simple for describing Huang et al.’s experimental data [4]. We conclude by showing that the minimal and original models make qualitatively similar predictions about the efficacy and robustness of treatment protocols not considered in our experimental dataset, which serves to further validate the sufficiency of the minimal model.

#### 2. Methods

We begin by introducing the previously developed mathematical model that describes tumor-immune dynamics subject to treatment with either OVs, DCs, or both [17, 18]. We then expound upon the variety of techniques that we utilized to address whether the model is minimally structured, given the experimental data it was designed to fit.

##### 2.1. Original Mathematical Model

Our original model contained the following six ordinary differential equations:

When all parameters and time-varying terms are positive, this model captures the effects of tumor growth and response to treatment with Ad/4-1BBL/IL-12 and DCs. By allowing various parameters and time-varying terms to be identically zero, other treatment protocols tested in Huang et al. can also be described (e.g., adenovirus that only mediates the release of 4-1BBL) [4]. These other models are defined more explicitly after the description of the full model.

In Equation (1), uninfected tumor cell volume, , grows exponentially (at a rate ), and upon being infected by an oncolytic adenovirus, , converts to infected cancer cell volume, , at a density-dependent rate (), where is the total volume of cells (tumor cells and tumor-targeting T cells, ). The tumor-targeting T cells indiscriminately kill both uninfected and infected tumor cells, with the rate of killing depending on the amount of IL-12 and 4-1BBL production (modeled through in the term ).

In Equation (2), newly infected tumor cell volume is accumulated at a rate of . The infected cells, , are lysed by the virus or other mechanisms (at a rate of ), thus acting as a source term for the virus by releasing free virions into the tissue space (with *α* virions released on average per cell, seen in Equation (3)). We again see the tumor-targeting T cells in action, killing (modeled in the term ).

In Equation (3), treatment with OVs is represented with the time-dependent term , determined by drug delivery and dosing schedules of interest. Due to lysis of , virions are released on average per cell lysed, as discussed earlier. Virions are not impacted by the T cell population—any loss of virions is due to natural decay, .

Equations (4) and (5) describe how the population of T cells, , and naïve T cells, , change in time. The activation and recruitment of tumor-targeting T cells can happen in three ways: (1) through stimulation of the naïve T cell pool, which at rate can asymmetrically divide to give rise to tumor-targeting T cells, due to increased IL-12 (and modeled as proportional to , at a rate of , as infected cells are the ones to release IL-12); (2) through stimulation of cytotoxic T cells due to increased 4-1BBL (also modeled as proportional to , at a rate of ); and (3) through production of T cells due to the externally primed dendritic cells, , (at a rate of ). T cells and naïve T cells also experience natural death (at rates of and , respectively).

Equation (6) describes how the population of injected dendritic cells changes over time. The time-dependent term represents the treatment with DCs, determined by the drug delivery and dosing schedule of interest. Dendritic cells decay at a rate of .

The data from Huang et al. was measured as tumor volume versus time for a variety of treatment protocols, averaged over 8-9 mice per protocol [4]. The treatment protocols increased in complexity, and therefore, our full model was designed and validated hierarchically against each dataset as follows:(i)*Model 1*: *No OV or DC treatment*. Under no treatment, tumor data are modeled using Equation (1) only, with all parameters other than set to 0.(ii)*Model 2*: *Treatment with OVs that replicate and lyse, but do not mediate the release of cytokines and costimulatory molecules*. Tumor data under this treatment are modeled using Equations (1)–(3) with and set to 0.(iii)*Model 3*: *Treatment with OVs that lyse tumor cells and*(a)*mediate the release of 4-1BBL (Ad/4-1BBL)*. Tumor data under this treatment are modeled using Equations (1)–(4) with and parameters in Equation (4) set to 0.(b)*mediate the release of IL-12 (Ad/IL-12)*. Tumor data under this treatment are modeled using Equations (1)–(5) with and parameters in Equation (4) set to 0.(c)*mediate the release of both IL-12 and 4-1BBL (Ad/4-1BBL/IL-12)*. Tumor data under this treatment are modeled using Equations (1)–(5) with parameter in Equation (4) set to 0.(iv)*Model 4*: *Treatment with DCs only*. Tumor data under this treatment are modeled using only Equations (1), (4), and (6). Equation (4) has and parameters set to 0 in this case, and Equation (1) has set to 0.(v)*Model 5*: *Treatment with Ad/4-1BBL/IL-12 and DCs*. This is modeled using the entirety of the system in Equations (1)–(6).

The experimental data increase from simple to more complex, allowing us to fit the model parameters in a hierarchical fashion. For instance, using Model 1, we fit the tumor growth rate to the control data in which the tumors grew without treatment. The best-fit value of was then used in subsequent versions of the model. Using previously fit parameters allowed us to reduce the dimension of parameter space at each step of the model development and fitting process. More details on the hierarchical development of the model can be found in our previous studies [17–19].

##### 2.2. 95% Credible Intervals and Local Sensitivity Analysis

In our original works [17, 18], a single best-fit value was determined for each parameter in system (1)–(6). In our later work [19], we expanded our understanding of the best-fit parameter values by identifying the potential distribution for each parameter. We did this by creating one thousand bootstrap replicates [23] from each of our experimental datasets (control, Ad only, Ad/4-1BBL/IL-12, etc.,) [19]. Each bootstrap replicate was created by sampling the mice in the original experimental dataset with replacement. The best-fit parameter values were found for each bootstrap replicate, and these were used to estimate the posterior marginal distribution on each fit parameter. In other words, each distribution we approximate is the distribution of the sample average of a parameter value, for one thousand populations of size . The interval in which we can be 95% certain that the true value of the parameter is found, the 95% credible interval, was then calculated from these approximated distributions. This is done by excluding the values that fall in the extreme tails of the distribution (2.5% of values in each tail). In this work, we identify poorly constrained parameters (that may not be needed in a minimal model) as those with very large 95% credible intervals.

We have previously performed a local sensitivity analysis on the full system in Equations (1)–(6), along with the submodels that describe the simpler treatment protocols detailed above [19]. The parameters we focus on are the ones whose values could not be readily ascertained from experiments, as detailed in [18]. For each submodel, our local sensitivity analysis entails performing an exhaustive search about the best-fit parameters and identifying all parameter sets that give a fit within 10% of the optimal fit.

The optimal fit is defined as the set of parameters that minimize the goodness-of-fit metric [19]:where is the average experimental tumor volume at day , is the tumor volume at day predicted by the submodel under consideration, and is the variance in the experimental tumor volume at day . The fractional term in Equation (7) is a dimensionless measure of the error in which the sum of the square error is divided by the variance in the experimental data. In this way, we require better fits to the average volume when the variance is small, in accordance with the principle of maximum likelihood estimation [24]. However, because instrumentation error in volume measurements is independent of tumor volume, calipers are imprecise for smaller tumor sizes [25, 26]. For this reason, we weigh the dimensionless term in Equation (7) by the average tumor volume, as this does not artificially bias the algorithm to fit well at small tumor sizes at the expense of fitting well over a majority of the data points [19].

In this work, we use these data to identify parameters that the model is highly insensitive to (parameters that can vary widely without affecting the goodness of fit of the model). We then consider fixing each of those parameters or even restructuring the equations to remove such parameters all together.

##### 2.3. Global Sensitivity Analysis

We next take a more holistic view of parameter space by performing a global sensitivity analysis using the Sobol method [27]. To do this, we view parameter space as our input space, and the tumor volume as predicted by the corresponding DE submodel in Equations (1)–(6) as the output space. The Sobol method allows us to determine how much of the total output variance is due to each individual parameter. If the Sobol Index corresponding to a particular parameter (defined below) is small, the model dynamics are insensitive to varying that particular parameter.

To use the Sobol method, we first determine a realistic domain for each parameter in the model. Define the vector such that each element of the vector corresponds to a parameter in our model, specified in Equations (1)–(6):

The full set of choices for each element of the vector defines the 14 dimensional parameter space ; though in all of our analyses, we only consider a subset of this vector. We assume a uniform distribution for each parameter in the parameter space, where the minimum and maximum value for each parameter is determined from its 95% credible interval.

Here, we define the total output space as all possible values of , where for a particular choice of parameters , we definewhere and are elements of the solution to the original model (1)–(6), for the particular choice of parameters with time to , here 30 days.

In order to determine how the output varies with each parameter, we consider how much each parameter contributes to the total variance of the output space [27, 28], . First we calculate the total variance of over the full parameter space, . If possible, we would next calculate the variance of , given the true value of , , to determine how much variance is lost due to being fixed. However, we do not know the true value of . Instead, we calculate the first-order sensitivity index (Sobol Index) for each of the parameters, , by calculating the expected value of over all other parameters conditioned on each value of : . We then compute the variance of the conditional expectation: .

is then defined as the amount that the parameter contributes to the total variance, normalized by the total variance:where, as previously stated, the range of values considered for each parameter is determined by its associated 95% credible interval.

Sobol’s key insight involves using an ANOVA decomposition to calculate the above variances [27, 28]. This allows for practical numerical simulations that employ Monte Carlo simulations to calculate the Sobol Indices [27]. In this study, we utilized an available MATLAB toolbox, GSAT, to calculate all of the reported Sobol Indices using the FAST method and Sobol distributions for sampling. (https://www.mathworks.com/matlabcentral/fileexchange/40759-global-sensitivity-analysis-toolbox). When Sobol’s method determines that varying a parameter over its domain causes to change minimally, will be small. Because such a does not contribute much to the overall variance of the output, the model is considered insensitive to the value of that parameter. On the contrary, when varying a parameter over its domain causes to greatly vary, will be large. A model is sensitive to parameters corresponding to large values.

##### 2.4. Parameter Sloppiness Analysis

It is well established that the quantitative behavior of multiparametric biological models is much more sensitive to changes in certain combinations of parameters than to others, a phenomenon known as “sloppiness” [29]. Herein, we will use an established methodology (see [29, 30] for more details) to identify parametric combinations to which the model is highly sensitive (i.e., stiff directions) and combinations to which the model is highly insensitive (i.e., soft directions). This is done by performing an eigenvector decomposition of the Hessian matrix corresponding to our cost function in Equation (7). Large eigenvalues correspond to stiff eigendirections, whereas small eigenvalues correspond to the soft eigendirections [29, 30].

The foundational step of performing a sloppiness analysis is to compute the cost manifold [30]. We do this by identifying all points in parameter space that give a fit within 10% of the optimal (generated for the previously-described local sensitivity analysis). The data are then normalized in two ways. First, the value of each parameter is scaled by the optimal value of that parameter, hence the normalized axes we see in Figure 1. Second, for each parameter set , we define the cost function to be the deviation in the value of at that parameter set compared with the value of at the optimal parameter set. As we are only interested in parameters that give a fit within 10% of optimal, each point shown in parameter space in Figure 1 has a value in the range 0 and 0.1. We then find the best-fit quadratic surface describing using the polyfit package for MATLAB (https://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn). Once the best-fit quadratic surface has been obtained (an example of which is visualized in Figure 1(a)), its associated Hessian matrix can be computed. The eigenvalues and eigenvectors of this matrix allow us to identify soft and stiff parameter directions.