Article Text
PDF
Abstract
Background Immunotherapies, driven by immune-mediated antitumorigenicity, offer the potential for significant improvements to the treatment of multiple cancer types. Identifying therapeutic strategies that bolster antitumor immunity while limiting immune suppression is critical to selecting treatment combinations and schedules that offer durable therapeutic benefits. Combination oncolytic virus (OV) therapy, wherein complementary OVs are administered in succession, offer such promise, yet their translation from preclinical studies to clinical implementation is a major challenge. Overcoming this obstacle requires answering fundamental questions about how to effectively design and tailor schedules to provide the most benefit to patients.
Methods We developed a computational biology model of combined oncolytic vaccinia (an enhancer virus) and vesicular stomatitis virus (VSV) calibrated to and validated against multiple data sources. We then optimized protocols in a cohort of heterogeneous virtual individuals by leveraging this model and our previously established in silico clinical trial platform.
Results Enhancer multiplicity was shown to have little to no impact on the average response to therapy. However, the duration of the VSV injection lag was found to be determinant for survival outcomes. Importantly, through treatment individualization, we found that optimal combination schedules are closely linked to tumor aggressivity. We predicted that patients with aggressively growing tumors required a single enhancer followed by a VSV injection 1 day later, whereas a small subset of patients with the slowest growing tumors needed multiple enhancers followed by a longer VSV delay of 15 days, suggesting that intrinsic tumor growth rates could inform the segregation of patients into clinical trials and ultimately determine patient survival. These results were validated in entirely new cohorts of virtual individuals with aggressive or non-aggressive subtypes.
Conclusions Based on our results, improved therapeutic schedules for combinations with enhancer OVs can be studied and implemented. Our results further underline the impact of interdisciplinary approaches to preclinical planning and the importance of computational approaches to drug discovery and development.
- computational biology
- oncolytic virotherapy
- clinical trials as topic
- drug therapy
- combination
- drug evaluation
- preclinical
Data availability statement
All Matlab code used for generating the virtual cohort and analysis of the cohort under enhancer injection multiplicity and VSV lag is freely available at the GitHub repository https://github.com/adriannejnner/In-silico-trials-for-combination-strategies-for-enhancing-vesicular-stomatisis-oncolytic-virus. The virtual cohort used for our analyses is also available in the repository. Data are available on reasonable request from the corresponding author (morgan.craig@umontreal.ca).
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
- computational biology
- oncolytic virotherapy
- clinical trials as topic
- drug therapy
- combination
- drug evaluation
- preclinical
Background
Oncolytic viruses (OVs) are genetically modified viruses designed to specifically target tumor cells.1 The antitumor effects associated with oncolytic virotherapy are mediated significantly by immune mechanisms, which can be either advantageous or disadvantageous depending on the type of virus.1 Although immunosuppression may improve viral oncolysis, this gain is achieved at the cost of antitumor immunity, a key factor for improving cancer therapies. The importance of considering immune-virus interactions is supported by the mechanism of action of the OV talimogene laherparepvec (T-VEC).2 T-VEC was the first US Food and Drug Administration (FDA)-approved OV and is a genetically modified form of herpes simplex virus that encodes the immunostimulatory cytokine granulocyte-monocyte colony-stimulating factor (GM-CSF). The OV’s effectiveness is amplified by its immunostimulatory counterpart,3 attesting to the need to find a reasonable balance between a multitude of immune mechanisms (such as viral clearance and the antitumor immune response) to achieve ultimate treatment success.
Therapeutic cancer vaccines are administered to cancer patients with the goal of eradicating tumor cells through strengthening the patient’s own immune response.4 Combination OV protocols or vaccination schedules use a sequential combination of immunologically distinct viruses to induce immunity, circumvent or mitigate the antiviral immune response, and ultimately enhance antitumor efficacy.5–7 Currently, there are three clinical trials investigating the efficacy of combining adenovirus and the OV Maraba as an anticancer vaccination treatment,8 9 two of which are in the USA (NCT02285816; NCT02879760) and one in Canada (control number 195876, protocol number: AD/MG1-MAGEA3-001).
The principal idea behind combination OV regimens is to stimulate alternate mechanisms of antitumor immunity that act cooperatively or synergistically to enhance therapeutic effect. In this way, OVs can be used in enhancer virus/primary virus regimens whereby the pre-existing immune response induced by the enhancer OV will improve the efficacy of the second primary heterologous OV administration.10 It has been shown that the development of an acquired antiviral immune response usually takes less than a week in treatment-naïve animals, leaving a small window of opportunity for oncolytic vectors to function.11 Accepting that the ensuing immune response dictates that viral oncolysis will inevitably be transient in nature, the anti-OV immune response can be usefully reoriented to enhance the therapeutic impact of the vector.
Bridle et al7 were among the first to investigate the synergy between combination OV therapies. They showed that the antitumor response to vesicular stomatitis virus (VSV) was weaker than the anti-VSV response. This led them to complement the initial injection of VSV with an injection of a different virus with the goal of harnessing the original anti-VSV response and improving the antitumor immune response. Bridle et al12 therefore, investigated VSV as a boost to adenovirus antigen. They found that VSV antigen produced a more tumor-specific CD8+ T cell response which was more cytotoxic in combination with adenovirus-antigen with increased cytokine and granzyme production.
The exact immune mechanisms through which OVs induce antitumor responses depend on the type of virus used and the transgenes encoded. Ilett et al1 showed that reovirus induced the priming of a CD8+, Th1-type antitumor response whereas VSV expression promoted a potent antitumor CD4 +Th17 response, and that priming with reovirus, followed by VSV significantly improved survival of B16 melanoma tumor-bearing mice versus virus alone. Previous work has also suggested that three low doses of adenovirus, followed by three low doses of vaccinia virus (VV) resulted in a superior antitumor efficacy versus six doses of either virus.13
Individually, VV and VSV have both been extensively investigated as possible oncolytic virotherapy agents.11 13–18 Morphologically and immunologically, VSV and VV are very distinct. VV is a complex double stranded DNA virus encoding a large number of genes with immune evading properties that allow the virus to establish local pockets of infection within an infected host at a tissue level.19 At the systemic level, VV is a highly immunogenic virus, eliciting strong T-cell-mediated and antibody responses.20 Due to the role VV played in the worldwide smallpox eradication program, it has long been recognized as an efficient therapeutic vaccine and has the longest and most extensive history of human use of any virus, which demonstrates its safety.20 In contrast, VSV is a genetically simple RNA virus (with only five gene products) that rapidly replicates and spreads within tumors. VSV is extremely sensitive to the antiviral effects of type I interferons (IFNs),19 which act to inhibit viral replication and spread in immunocompetent (IC) hosts.21 22
The feasibility of using VV and VSV together in combination OV treatment was previously demonstrated, with the potential to improve therapeutic outcomes in triple negative breast cancer (TNBC).19 Le Boeuf et al19 used VV naturally expressing the viral gene product B18R, an IFN receptor decoy that locally antagonizes the cellular antiviral response initiated by type I IFNs, in parallel with a recombinant version of VSV expressing fusion-associated small-transmembrane protein to further enhance VV’s ability to spread through an infected monolayer. The combination of these viruses resulted in a ‘ping pong’ oncolytic effect wherein VV enhanced the ability of VSV to replicate and/or spread in tumor cells. In their work, Le Boeuf et al19 only considered a single administration of the combined dosage protocol (VV+VSV). A rational approach leveraging quantitative, predictive modeling and experimental results would help to delineate the therapeutic potential of combined enhancer VV with VSV, and further preclinical investigations into combined oncolytic virotherapy strategies.
The translation of OVs from preclinical studies to clinical implementation is a major challenge. Solving this obstacle requires answering fundamental questions about how to effectively design and tailor schedules to provide the most benefit to patients. Here, mathematical and computational biology help to identify strategies that offer durable therapeutic benefits prior to human trials.23–32 Interpatient heterogeneity is a defining obstacle in cancer therapy, and patient-to-patient variability in cancer can cause finely tailored treatment protocols to exhibit extreme disparate antitumor responses.33 34 Quantitative approaches have similarly been leveraged to integrate experimental data and identify robust optimal treatment protocols,33 with quantitative systems pharmacology models contributing to decision making at the regulatory level.35
In particular, virtual clinical trials36 37 (or in silico ‘twins’38 39) have recently been used in preclinical research to make ‘go or no go’ decisions.33 40–46 We have previously developed a computational biology model describing tumor-immune interactions and systemic cytokine concentrations over time,47 which we used to determine the optimal combination of GM-CSF and OV.46 We predicted that appropriately eliciting immune responses could significantly improve 5-year patient outcomes. Jafarnejad et al48 conducted an in silico clinical trial of anti-PD-1 molecule nivolumab for non-small-cell lung cancer calibrated to human patient clinical trial data. They predicted that patients with adjuvant nivolumab treatment in addition to the clinical trial protocol of neoadjuvant nivolumab treatment, followed by resection produced a durable response. With a focus on cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), programmed cell death protein 1 (PD-1) and programmed death-ligand 1 (PD-L1) blockade in melanoma, Milberg et al49 similarly leveraged a virtual clinical trial to predict the performance of therapeutic combinations given heterogeneous patient characteristics. Using their model validated to measurements from clinical trials, they predicted that response rates were higher for anti-PD-1/PD-L1 vs anti-CTLA-4/PD-1 combinations, and that anti-PD-1 administered before anti-CTLA-4 produced a greater response than the converse, consistent with clinical results. Applications of virtual patient ‘twins’ are not only specific to oncology but have also been applied in drug development, for example, to estimate the long-term effects of a treatment from short-term placebo-controlled trial measurements.36 39
Here, we expanded our mathematical model to investigate the therapeutic potential of enhancing VSV efficacy using VV as an enhancer. Using an in silico trial platform, we leveraged our computational biology model to predict tumor burden under clinically actionable combination OV-therapy administration schedules. Motivated by the results of Le Boeuf et al19 where VV was shown to enhance the efficacy of VSV, we then interrogated the impact of enhancer (VV) multiplicity and the lag in VSV administration on tumor growth dynamics to establish an optimal enhancer number and VSV lag. Our results suggest that intrinsic tumor characteristics, mainly tumor aggressivity, are the primary drivers of therapeutic response and ultimate success. Importantly, we show that these attributes can be exploited for patient stratification and to tailor therapeutic protocols.
Methods
Mathematical model of combination OV therapy and response
We extended our previous model for tumor growth and resistance to treatment with virotherapy and the resulting immune response46 to consider the impact of using dual VV ‘enhancer’ and VSV OV injection. In this scenario, the antitumor immune response is either upregulated or downregulated depending on the type of virus injected. We considered total vaccinia and vesicular stomatitis virions (
and 
, respectively) and two corresponding infected cell populations Extra close brace or missing open brace and Missing \left or extra \right. Parameters from Cassidy and Craig46 relating to viral kinetics were taken here to be virus dependent (subscript notation for 
and ω). Additionally, we included an immunomodulation term 
that modulates the production of immunostimulatory cytokines. The complete set of model equation is provided in the online supplemental technical information. A summary of the biological assumptions and model schematic is given in figure 1.
Supplemental material
Schematic representation of the tumor growth model under combination OV-therapy. (A) Biological assumptions for the combination OV-therapy interactions between VV (enhancer) and VSV oncolytic viruses. Infection of cells by either VV or VSV results in cell lysis, whereby new virus particles are released along with a cocktail of antigens, antivirals and cytokines. For simplicity, we considered each virus to release associated cytokines concentrations that can independently instigate the recruitment of immune cells (such as phagocytes). However, in the presence of both VV and VSV infections, we assume the cytokine production decreases the recruitment of immune cells, allowing for a more targeted immune response and virus-induced cell lysis. Additionally, VSV releases antivirals that block the intracellular replication of the virus and the infection of neighboring cells. In comparison, lysis of VV infected cells produces antivirals that downregulate the antivirals produced by VSV, allowing for infection and replication to occur. Once activated, these immune cells induce cell apoptosis of uninfected cancer cells. (B) In the model, quiescent susceptible cells (light blue) activate and begin division by transitioning into the 
phase of the cell cycle. Cells exit 
to enter the active phase (mitosis) and complete division. Most susceptible cells in the active phase re-enter quiescence after mitosis, however, certain dividing cells may mutate into an immune-resistant lineage (red). Immune interactions are driven by immune cells who encounter quiescent
and actively dividing susceptible tumor cells. Tumor-immune interactions increase proinflammatory cytokine concentrations to recruit additional immune cells to the tumor site. VSV and VV infect both normal and immune-resistant tumor cells, creating virus-specific infected cell pools. These infected cells undergo lysis releasing new virus progeny. The virus also influences the cytokine production which controls the immune cell production and activity. VSV: vesicular stomatitis virus; VV: vaccinia virus.
Briefly, the following biological interactions were added to the model described in Cassidy and Craig46:
VSV and VV are morphologically and cytotoxically distinct19 and therefore have the following virus-specific characteristics: virion-cell infection rates (
and ), virion induced cell lysis rates ( and ), virion burst sizes ( and ), and virion death rates ( and ).Viruses modulate the production of cytokines through either promoting an inflammatory or anti-inflammatory immune response,50 modeled through an immune modulation constant (
) that controls the rate of cytokine production from the immune interaction with cycling tumor cells (infected and uninfected). As , cytokine production is reduced, recapitulating an anti-inflammatory immune regulation, and as , normal inflammatory immune response is recovered.VV downregulates the production of antiviral factors, which aids the spread of VSV.19 Therefore, we consider VV to downregulate cytokine production, instigating an anti-inflammatory response (ie,
), whereas VSV upregulates cytokine production from both infected and uninfected tumor cells (ie, ).Limitations on the binding of immune cells to cognate growth factors or signals50 due to the simultaneous infection of both VV and VSV in the tumor were represented by higher production of inflammatory cytokines (
) and a lower rate of maximal immune cell production () when VSV particles are introduced after VV has commenced replication.
All other interactions are as in Cassidy and Craig.46 For more information on the biological interactions and their model implementation see online supplemental technical information.
Estimation of vaccinia (VV) and VSV viral and immune related parameters
Parameters of the model were estimated via a hierarchical fitting algorithm in which subsets of the model were fit to different experiments using VV and VSV. Full details are provided in the online supplemental technical information. Briefly, we sequentially fit model parameters to the experimental measurements from HT29 and 4TI cell lines from Le Boeuf et al19 and Rausch et al51 (online supplemental figures TS1–TS3). Tumor growth parameters were obtained by fitting the rate quiescent cells enter the 
phase (
), the rate HT29 (human colorectal adenocarcinoma) and 4T1 (murine mammary carcinoma) cells leave 
to enter the active phase (
), and the rate cells undergo apoptosis in 
phase (
), in immunodeficient (ID) mice (ie, HT29/ID and 4T1/ID obtained by Le Boeuf et al19 and Rausch et al,51 respectively). We then fixed these parameter values and estimated the viral kinetic parameters 
and 
from VSV, VV and VV+VSV treated HT29 tumor growth measurements in ID mice (ie, HT29/ID-VSV, HT29/ID-VV and HT29/ID-VV+VSV obtained by Le Boeuf et al19). Then, using IC mice experiments from Rausch et al,51 we estimated the immune cell-tumor cell contact rate 
, the immune cell digestion constant 
, the cytokine production half effect 
and the maximal immune cell production rate 
. To account for the effects of humoral immune responses the viral-kinetic parameters 
were then recalibrated to the presence of the immune system using 4T1 tumor growth measurements in IC mice under VSV, VV and VV+VSV (ie, 4T1/IC-VSV, 4T1/IC-VV and 4T1/IC-VV+VSV obtained by Le Boeuf et al.19).
Supplemental material
Generation of in silico individuals and patient cohorts
To reflect interindividual variability and the heterogeneity in treatment outcomes, we generated a unique set of parameters to represent individual virtual patients (figure 2A, no human patients were involved in this study). For this, we sampled tumor and immune cell-related parameters 
(where τ is the expected tumor cell cycle duration) from a normal distribution with mean 
corresponding to the parameter value returned in the hierarchical fitting described above.
In silico trial strategy recapitulates experimentally observed variability. (A) (1) Model parameters are established by calibrating experimental results to the model’s predictions. A distribution of responses centered at the mean of the experimental data is then used to generate parameter sets representing virtual individuals. (2) To populate the trial, each virtual patient’s tumor growth is simulated to determine whether they are candidates for the trial. Patients whose tumor growth is acceptable (ie, clinically relevant) are placed into repeated identical cohorts. (3) Alternative treatment schedules are then tested on each cohort by simulating individual virtual patient responses with the mathematical model and summarizing cohort level outcomes (such as mean and SD of responses). (4) optimal actionable schedules are then inferred by comparing cohort level and individual outcomes. (B) Tumor growth (relative to tumor volume on day 6) over time in absence of treatment. Black line: model fit; red stars: experimental observations measured by Le Boeuf et al,19 gray shaded region: distribution of growth from full cohort of patients. (C) Virtual patients were ordered based on intrinsic tumor growth rates r (top and bottom 10% denoted by shaded regions). OV: oncolytic virus.
To avoid the inclusion of non-realistic virtual individuals, we verified that each parameter set resulted in a tumor growth within two SD of the experimental measurements and the mean prediction at each corresponding data time point (figure 2B). Using this approach, we created 200 patients with parameter values normally distributed about the mean empirical or fitted value (online supplemental figure S1, online supplemental information), rejecting 265 parameter sets for not meeting the inclusion criteria. Since the maximal immune cell production rate (
) changes when VSV is introduced, we assumed each individual patient’s parameter varied equivalently. To recapitulate heterogeneity in initial tumor size and immune populations, we simulated an initial seeding of 
tumor cells, along with an initial cytokine concentration 
and immune cell count 
(day 0) for each in silico patient parameter set and fixed the initial conditions for treatment to be the tumor and immune populations on day 6 (online supplemental figure S2).
Kaplan-Meier survival curves for each cohort treated using the combination protocols were used to compare the effectiveness of the different trials. To determine a cull threshold for the survival of the virtual patients, we extrapolated the Kaplan-Meier curves in Le Boeuf et al19 by taking populations of 10 randomly sampled individuals and calibrating to their cumulative survival curve, giving a volume threshold (online supplemental technical information, online supplemental figure TS4). A local parameter sensitivity analysis of the mean empirical values showed that tumor growth was closely related to cell cycling rates and immune stimulation (online supplemental figure S3, online supplemental information).
We also quantified the growth rate r of the control tumors by approximating the growth curves with an exponential growth function from day 12 to 18, that is, 
, to obtain an estimate of later tumor growth. This measurement period was chosen to account for the discernable differences between cell lines after day 12 and the experimental end point at day 18.19 The implicit parameter r describes the aggressivity of the tumors, with high r corresponding to aggressively growing tumors and low r corresponding to slowly growing tumors. Ordering patients by their tumor growth rate showed a gradual increase in r values across the cohort (figure 2C). Numerical simulations, the creation of the virtual cohort, all statistical analyses and figure generation were performed using Matlab R2019b.
Results
Tumor aggressivity dictates the optimal number of enhancer injections
The effect of the multiplicity of VV enhancer injections on the success of therapy is largely unknown, especially for a cohort with varying underlying tumor growth rates. To assess the impact of the number of enhancer administrations, we simulated our virtual cohort after multiple daily enhancer infections followed by a VSV injection given 7 days after the last enhancer (figure 3A). The total number of VV injections (
) ranged from 1 to 7 and the dosage sizes were fixed to those used in the Le Boeuf et al19 (online supplemental technical information).
Influence of enhancer injection multiplicity on tumor burden. (A) The effects of enhancer multiplicity (
) were investigated by simulating 1–7 VV enhancers, with VSV administered 7 days after the final enhancing dose. Tumor growth was assessed 15 days after the administration of VSV. (B) Distribution in number of tumor cells 15 days after VSV administration with respect to the multiplicity of enhancers. Central mark (red) indicates median, bottom edge denotes the bottom quartile, top edge denotes the top quartile. Significance indicators report the non-significant results of a Kolmogorov-Smirnov test for significance of difference between distributions 
. (C) Tumor growth dynamics from last enhancer to 15 days after VSV administration for protocols with 1 (blue) and 7 (red) enhancers. Mean is denoted by a solid line, SD by shaded regions of same color and individual virtual patient values are plotted as circles. (D) Kaplan-Meier survival curves for protocols with 1 (dark blue), 2 (light blue) and 7 (red) enhancers. No significant difference between protocols from 2 to 7 enhancers was found (measured by log-rank test for significance, 
). VSV: vesicular stomatitis virus; VV: vaccinia virus.
To assess the performance of different enhancer schedules, we quantified the number of tumor cells 15 days after VSV administration (VSV +15 days; figure 3B). In particular, we found that the number of VV enhancers (
) did not significantly impact the average tumor size at VSV+15, and that the distribution of tumor sizes between sequential enhancer multiplicities, that is, one enhancer and two enhancers, two enhancers and three enhancers etc., was not significantly different. There was, however, a difference in the distribution of tumor sizes between one enhancer and all other enhancer multiplicities from three enhancers onwards (as confirmed by a Kolmogrov-Smirnov test, 
, see online supplemental figure S4A for significant pairings). We also determined that the mean number of tumor cells was only significantly different between one and six enhancers, and one and seven enhancers (pairwise t-test, 
and 
, respectively; online supplemental figure S4B, online supplemental information). This suggests that the duration of treatment and tumor aggressivity could be the principle drivers of the distribution in tumor sizes at 20 days, irrespective of the number of enhancers.
Nonetheless, we were able to distinguish low and high responders by their tumor growth at VSV +15 days after 7 vs 1 enhancer (figure 3C; results for all multiplicities of enhancer injections in online supplemental figure S5A). Though the mean dynamics of the number of tumor cells are qualitatively similar after one or seven enhancers, there is a statistically significant difference in the mean number of tumor cells at VSV +15 days for one enhancer compared with seven enhancers. The difference in variance of cohort responses can be explained by the fact that the last enhancer is administered on day 1 for the one enhancer protocol vs day 6 under the seven enhancer schedule, implying that tumors have ultimately been growing for a longer period of time under the seven enhancer protocol, supporting the conclusion that therapeutic success is largely driven by intrinsic aggressivity (ie, higher growth rates). Log-rank tests of the cohort’s Kaplan-Meier survival curves also showed that the one enhancer protocol was significantly different from all others (figure 3D and online supplemental figure S5B in the online supplemental information).
An increase in the variability of responses was observed with a corresponding increase in the number of enhancers (figure 3B). To investigate whether there was a link between a patient’s intrinsic tumor growth rate r and the optimal number of enhancers, we established each patient’s optimal number of enhancers through numerically simulating all possible treatment protocols and finding the one that minimized tumor burden at VSV +15. We found that the optimal number of enhancers grew with decreases in intrinsic tumor growth rates (figure 4A). Further, our results show that individuals with high growth rates consistently had worse outcomes, even during ‘optimal’ treatment with two enhancers. To discern patterns of responses in less aggressive tumors and higher growth rates, we ordered the number of enhancers from best protocol to worst protocol for each patient, based on the tumor size 15 days after the VSV injection (figure 4B). Patients with low tumor growth rates were found to perform best under treatment with seven enhancers and worst under treatment with a single enhancer, whereas patients with high tumor growth rates performed best with two enhancers and worst under seven enhancers. These results suggest that there is a significant difference in the efficacy of protocols given a patient’s intrinsic tumor growth rate.
Individual responses to multiple enhancer injections protocols are stratified by intrinsic tumor growth rates. (A) Waterfall plot of the change in tumor size 15 days after VSV administration. Bar color depicts the optimal number of enhancers. Corresponding tumor growth rate (r) for each patient are plotted in gray. (B) Ordering of protocols from best (bottom row) to worst (top row) for each patient based on the tumor size 15 days after VSV administration. Corresponding tumor growth rates are plotted above (patient ordering identical based on intrinsic tumor growth rate as in A). VSV: vesicular stomatitis virus; VV: vaccinia virus.
Shorter VSV lags are necessary to slow aggressive tumors
Thus far, enhancer multiplicity was investigated with VSV administered 7 days after the last enhancer. However, the number of days between the final enhancer and the VSV administration can also influence the priming of the immune response and the efficacy of the treatment. To measure the impact of the VSV lag (
days), we simulated a single administration of VSV 1–15 days after the final enhancer, with dosages fixed as in Le Boeuf et al19 (figure 5A). The number of enhancers was set to either 1 or 7, based on the enhancer multiplicity results described previously. As before, therapeutic response was assessed by the number of tumor cells (ie, tumor size) at VSV +15 days (figure 5B).
Longer VSV lags have a detrimental effect on therapeutic success. (A) Inspired by the results for enhancer multiplicity, the effects of the length of VSV lags were investigated by simulating either a single enhancer protocol (left) or a seven enhancer protocol (right). VSV lags were varied from 1 to 15 days after the final enhancer. Tumor growth was assessed 15 days after the administration of the last VSV. (B) Distribution in number of tumor cells 15 days after VSV administration with respect to the duration of VSV lag after 1 (orange) or 7 (blue) enhancers. Central mark (red) indicates median, bottom edge denotes the bottom quartile, top edge denotes the top quartile and individual virtual patient values are plotted as circles. (C) Tumor growth as a function of time. Markers indicate significant differences in final tumor size (t-test, 
). (D) Cytokine concentrations as a function of time. (E) Immune cells as a function of time. In C, D, E: 1 enhancer, 1-day VSV lag (solid blue), 7 enhancers, 1-day VSV lag (dashed blue), 1 enhancer, 15-day VSV lag (solid red), 7 enhancers, 15-day VSV lag (dashed red). VSV: vesicular stomatitis virus; VV: vaccinia virus.
Similar to what was observed with the enhancer multiplicity, increasing the VSV lag increased the variance of the cohort’s response (figure 5B). A VSV lag of 1 day, irrespective of the number of enhancers, produced the smallest average tumor size, with a very small distribution of responses at VSV +15 days. In comparison, we observed more dispersion in overall responses for VSV lags of 7 days or more, with some individuals achieving lower tumor sizes than for shorter VSV lags (see online supplemental information, online supplemental figure S6A).
Over longer time, the separation in the tumor size under a 1-day vs 15-day VSV lag becomes clearer (figure 5C). On average, 15-day VSV injection-lags performed the worst of all tested scenarios, especially at longer time points (65 days past VSV). This further solidifies that, as opposed to enhancer multiplicity, the time to VSV administration is the key determinant of tumor size. Since the immune cell population is essentially stabilized by VSV +15 days (figure 5E), and cytokine concentrations are saturated (figure 5D), the dynamics that occur immediately after the VSV injection must lead to the divergence of long-term tumor behavior.
As we hypothesized that interindividual variability would significantly impact treatment responses, we next investigated the optimal VSV lag for each individual after either one or seven enhancers. For the one enhancer protocol, the optimal VSV lag was a single day for all but three individuals with particularly slow growing tumors, for whom a 5-day VSV lag was best (figure 6A). However, we expect this response to be not particularly significant, as the difference in tumor cell numbers between the optimal and the second most optimal protocol for these three individuals was negligible compared with the rest of the cohort (results not shown).
Optimal VSV lag is stratified by intrinsic tumor growth rates. (A) Waterfall plot of the change in tumor size 15 days after the last VSV administration for the one enhancer protocol. Bar color depicts the optimal VSV lag. Corresponding tumor growth rate (r) for each patient are plotted in gray. (B) Waterfall plot of the change in tumor size 15 days after the last VSV administration for the seven enhancer protocol. Bar color depicts the optimal VSV lag. Corresponding tumor growth rate (r) for each patient are plotted in gray. (C) Ordering of protocols from best (bottom row) to worst (top row) for each patient based on the tumor size 15 days after the last VSV administration for seven enhancer protocol. Corresponding tumor growth rates are plotted above (patient ordering identical based on intrinsic tumor growth rate as in B. The ordering of the optimal VSV lag for the one enhancer protocol is provided in online supplemental figure S6 in online supplemental information). VSV: vesicular stomatitis virus.
Tumor responses to treatment were increasingly stratified for optimal schedules after seven enhancer administrations (figure 6B), likely related to the increased interindividual variability observed in this case. We found that optimal schedules for virtual patients with the slowest intrinsic growth rates required a 15-day VSV lag vs 1 day for those with slow growing tumors. Indeed, there was a clear delineation between aggressively growing tumors that required as short of a lag as possible for maximal therapeutic effect and slower growing tumors that responded best with as long of a VSV lag as possible (figure 6C and online supplemental figure S7 in the online supplemental information). Crucially, individuals with the slowest tumor growth were predicted to have the most meaningful responses, completely recovering in some cases. This again underlines that patient stratification and schedule tailoring is crucial for ensuring the most meaningful clinical response.
Individualizing enhancer-VSV scheduling
To further delineate vaccination scheduling in distinct subcohorts, we tailored VV enhancer multiplicity (between 1 and 7) and VSV lag (between 1 and 15 days) for each virtual individual according to intrinsic tumor characteristics (figure 7A). The individual optimal protocol was determined by simulating all scheduling possibilities and minimizing tumor size at VSV +15 days.
Individualized schedules are determined by tumor aggressivity and risk stratification according to tumor aggressivity is necessary for optimal outcomes. (A) Optimal number of enhancers (yellow), optimal VSV lag (fuchsia) and relative tumor size 15 days after last VSV versus untreated control (purple) as a function of intrinsic tumor growth rate. For all but a subset of the least aggressive tumors, individualized protocols called for a VSV lag of 1 day, with fewer than seven enhancers. (B) Two new cohorts of patients were generated with either aggressive tumor growth (
, purple) or slow tumor growth (
, green). Original Le Boeuf et al data (red stars) and model fit to the original data (black curve) as in figure 2C. (C) Each cohort was simulated according to the previously determined optimal aggressive protocol (one enhancer followed by a VSV 1 day later) and optimal slow protocol (seven enhancers followed by a VSV 15 days later). To assess the effect specificity of each protocol, a cross-over trial wherein virtual patients with fast growth were treated with the slow protocol and vice versa was performed. (D) Kaplan-Meier survival curves for the two cohorts under the two different protocols. (E) To confirm the robustness of the aggressive and slow protocols, the optimal number of enhancers and VSV lag was then determined for patients in the new cohorts. The results of the newly generated cohort were then compared with the original cohort in A. Tumor size 15 days after the VSV was assessed. Original cohort individualized therapy compared with the new slow growth (top left) and new aggressive growth (top right) individualized schedules. Overlays of the corresponding optimal number of enhancers and VSV lag for each patient from old protocol versus the new slow growth (bottom left) and aggressive growth (bottom right). VSV: vesicular stomatitis virus; VV: vaccinia virus.
As before, we found clear stratification between optimal protocols for slowly growing (seven enhancers followed by a 14 or 15 days VSV lag) and fast growing (one enhancer followed by a 1-day VSV lag) tumors (figure 7A). Interestingly, for a range of low intrinsic growth rates r, optimal schedules resulted in near complete tumor removal, whereas we found a jump in the number of tumor cells, followed by a linear dependence on the intrinsic growth rate after a critical value around 
1/day. Clinically, an r value of 0.03 corresponds to a tumor doubling time of 23 days, which was above the original cohort average of 15 days (figure 2C). Overall, these results underline that tumor aggressivity is the determining factor for combination OV scheduling and the outcome of enhancer-VSV therapy.
Cross-validating tailored strategies
To confirm that outcomes are improved by employing therapeutic strategies based on tumor aggressivity and investigate the robustness of our stratification strategy, we generated two new cohorts comprised of 50 virtual individuals with slow growing tumors (
) and 50 with aggressively growing tumors (
). These ranges for r where chosen to correspond to the top and bottom 
of the initial cohort’s r value (figures 2C and 7B). We next simulated the tumor growth of these cohorts under the previously determined optimal protocols, that is, 7 enhancers and a 15-day VSV lag for slow growing tumors, and one enhancer and a 1-day VSV lag for aggressive tumors (figure 7C). Additionally, we simulated each cohort under the alternate optimal protocol, that is, aggressively growing tumors with the slow tumor growth protocol and vice versa.
Overall, in terms of survival and irrespective of protocol, the aggressive tumor growth cohort performed markedly worse than the slow tumor growth cohort (figure 7D). As these represent the 10% most aggressive tumors of the original cohort, it is not surprising that the efficacy of therapy is minimal. In contrast, survival in the slow growing cohort was markedly different under the two protocols: when treated with their optimized protocol, all individuals survived, and while survival declined when treated with the aggressive tumor growth protocol, it remained overall stronger than in the aggressive tumor cohort treated with its matched optimal treatment strategy. Nonetheless, both strategies perform better than the no treatment case (results not shown). To further confirm the optimality of the aggressive and slow tumor growth protocols, we also determined each virtual patient’s optimized combination schedule for the new cohorts (figure 7E). Unsurprisingly, the optimal protocols were the same as in the larger original cohort, implying there is a robust link between tumor aggressivity and the optimal combination OV-therapy protocol. Interestingly, for patients with extremely slow growing tumors (r close to 0.0196, doubling time of roughly 35 days) the optimal VSV lag was slightly shorter (
days) than the rest of the slow growing tumor cohort. Thus, we posit that a combined OV therapeutic vaccination approach with VV +VSV will be most effective for slow growing tumors.
Discussion
In 2015, the US FDA approved T-VEC for the treatment of non-resectable late-stage melanoma, making it the first OV to reach the Western market. However, despite much promise, the efficacy of OV monotherapy has been limited.52 53 In response, combined OV schedules hold much promise as an effective cancer therapy capable of eradicating tumor cells through virus infection, immune recruitment, and by providing a long-lasting durable response. Results from combined OV strategies are encouraging, with three clinical trials underway for an adenovirus and OV Maraba anticancer combination OV treatment.8 9
A major obstacle to the clinical implementation of combination OV-therapy protocols is designing promising and optimal therapeutic schedules. Further, the reproducibility of protocol efficacy must be demonstrated in heterogeneous patient cohorts. For this preclinical planning, mathematical and computational biology have a large role to play in predicting therapeutic responses and designing effective strategies. Leveraging our previous computational model,46 we developed an in silico model of combination OV-based therapeutic vaccination with vaccinia (VV) and VSV OVs to test the heterogeneous response to and optimality of an enhancer virus and VSV protocols. Each generated virtual patient was created based on a realistic distribution of model parameters, with growth following the trend of experimental results.
We found that the number of enhancers does not significant impact the average response of our generated virtual cohort. Though perhaps unintuitive, this is likely due to a saturation in the initial immune response. Investigating this further, our results show that while the variance of tumor sizes increased with the number of enhancers, the overall survival of the cohort did not vary significantly. Ultimately, no single optimal protocol was found for most of the cohort. However, at the individual level, there was a significant difference in outcomes found when optimizing the number of enhancers: for tumors with low intrinsic growth rates, a larger number of enhancers is necessary to be effective, whereas aggressive tumors required fewer enhancers. The latter finding supports the idea of ‘hitting hard, hitting early’ for fast growing tumors. This clear stratification based on tumor aggressivity suggests that the effectiveness of the enhancer-VSV protocol is largely a function of the relationship between viral replication and tumor growth. For a given initial tumor size, faster growing tumors will have more cells for the virus to infect and subsequently lyse, so there is a trade-off between the number of enhancers and the delay in administering the VSV to ensure that cells are sufficiently infected and subsequently recruiting immune cell subsets. On the other hand, for slowly growing tumors, more enhancers are needed to load the tumor microenvironment with virus and sufficiently activate the immune system. An advantage of OVs compared with other immunostimulatory compounds is that the concentration of OVs in the tumor microenvironment will initially increase due to viral replication, whereas other drugs will experience rapid clearance.
Conceptually, longer VSV lags should increase the mean tumor size, given that tumors are not eradicated by the priming protocol and are thus continuing to grow prior to administration of VSV. Indeed, we predicted that the optimal VSV lag should overall be between 1 and 4 days, or shorter than the currently used 7 days. Similar to what was observed for enhancer multiplicity, we found that increasing the VSV lag increased the dispersion in outcomes, irrespective of the number of enhancers. Thus, it is not necessarily the number of enhancers or VSV lag driving the variation in treatment response, but rather the duration of the treatment.
While there are currently no clinical trials combining VV and VSV for the treatment of TNBC, there are other clinical and experimental results that support the findings of our virtual clinical trial. A handful of clinical trials have been conducted with VV or VSV individually,9 54–56 and there is currently a phase I trial in stage III-IV melanoma using VSV-IFNbetaTYRP1 (NCT03865212). The randomized phase 2 dose-finding trial of Pexa-Vec (an oncolytic VV with a gene encoded to increase expression of GM-CSF) for the treatment of on advanced hepatocellular carcinoma found a clear distinction in survival based on dosage, with the high-dose cohort having an increase in overall survival compared with the low-dose cohort.54 Further, the high-dose protocols performed equivalently on patients with or without metastases, whereas, the low dose protocol was only optimal for patients without metastases. As the presence of metastases is an indicator of disease aggressivity, these findings align with our conclusion that the stage and overall aggressiveness of a patient’s disease is a determinant for oncolytic VV treatment efficacy, and underlines the importance of patient stratification into appropriate schedules for oncolytic VV therapeutic success.
Further, high levels of Ki-67 expression in TNBC have been found to correlate strongly with more aggressive proliferation and poor prognoses.57 By evaluating a total of 1800 patients with early invasive TNBC, Zhu et al57 determined that adjuvant chemotherapy was associated with better overall survival in patients with higher Ki-67 expression than patients with low Ki-67 expression, with adjuvant chemotherapy having no effect on disease-free survival in the latter group. Given that Ki-67 is a marker of proliferation, these data corroborate our finding that tumor aggressivity can be a predictor of treatment success in TNBC, and that using Ki-67 expression as a threshold for therapeutic planning and prognostic factor may improve survival.
There are certain limitations to our model formalism. Specifically, the stimulated immune response only impacts the growth of the tumor and only has a secondary effect on the virus: reducing the number of cells available for the virus to infect and lyse. Despite this shortcoming, our model was able to replicate the observed dynamics of VV and VSV in IC mice and accommodate for the anti-inflammatory and proinflammatory responses to these viruses. Future iterations of the model could build on this and develop a more complex model of the immune response to combination OV-based vaccine therapy. In doing so, our virtual clinical trial platform could be used to optimize combined OV-immunotherapy (such an anti-PD-1 immune checkpoint inhibitor1) and investigate whether we see a similar segregation of the optimized protocol based on tumor aggressivity. In addition, future experiments investigating combined VV and VSV treatment in other tumor lines or humans will allow for further model validation.
Through rational considerations, we developed a quantitative approach to therapeutic cancer vaccination that provides actionable and clinically relevant scheduling recommendations that can be easily translated from bench to bedside using complementary methodologies. Current experimental work in therapeutic vaccinations could provide effective novel cancer therapeutics1 58 and there is a growing push towards personalized tumor-specific vaccination therapies.59 Unfortunately, the pace of immunotherapy innovation is limited by clinical trial requirements. Here, we put forward a strategy to fill that gap, helping to define the next phase of combined OV regimens.
Data availability statement
All Matlab code used for generating the virtual cohort and analysis of the cohort under enhancer injection multiplicity and VSV lag is freely available at the GitHub repository https://github.com/adriannejnner/In-silico-trials-for-combination-strategies-for-enhancing-vesicular-stomatisis-oncolytic-virus. The virtual cohort used for our analyses is also available in the repository. Data are available on reasonable request from the corresponding author (morgan.craig@umontreal.ca).
Ethics statements
Patient consent for publication
Acknowledgments
All authors thank Dr Jean-Simon Diallo for correspondence regarding the published experiments from Le Boeuf et al.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
Twitter @adrianne_jenner, @tcass20, @lab_craig
Correction notice This article has been corrected since it first published. The provenance and peer review statement has been included.
Contributors ALJ, TC, M-CB-D and MC conceived of the study. ALJ, KB, TC and MC constructed the model. ALJ performed simulations. ALJ and MC wrote the article. All authors approved the manuscript.
Funding ALJ: Fonds de recherche Santé Québec Programme de bourse de formation postdoctorale pour les citoyens d’autres pays, Centre for Applied Mathematics in Bioscience and Medicine (CAMBAM) Postdoctoral fellowship. TC: Natural Sciences and Engineering Research Council of Canada (NSERC) PGS-D, U.S. Department of Energy 89233218CNA000001, NIH grants R01-AI116868 and R01-OD011095. M-CB-D : Institut du cancer de Montréal, Fonds de recherche Québec Santé, Fondation du cancer du sein du Québec. MC: NSERC Discovery Grant and Discovery Launch Supplement RGPIN-2018-04546.
Competing interests None declared.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.