Standing Variations Modeling Captures Inter-Individual Heterogeneity in a Deterministic Model of Prostate Cancer Response to Combination Therapy

Simple Summary Studying rare outcomes in cancer is challenging because observation of the rare event may require a very high number of patients or experimental animals. Here, we propose a new, predictive approach to understanding the biological mechanisms underlying such rare events in cancer treatment outcomes. We take as a case-study, the treatment of metastatic, castration-resistant prostate cancer with the live cell vaccine, sipuleucel-t (Provenge). Only a fraction of patients benefited from Provenge; that is, clinical success is a rare event. It remains an open question why Provenge conferred such a modest survival benefit. Our modeling paradigm captures the inherent heterogeneity that characterizes individuals in a population, and provides an explanation for the observed clinical outcomes of treatment with Provenge. Our approach readily generalizes to a range of emerging cancer immunotherapies, and more generally, to predicting and understanding how a population responds to any intervention targeting a human disease. Abstract Sipuleucel-T (Provenge) is the first live cell vaccine approved for advanced, hormonally refractive prostate cancer. However, survival benefit is modest and the optimal combination or schedule of sipuleucel-T with androgen depletion remains unknown. We employ a nonlinear dynamical systems approach to modeling the response of hormonally refractive prostate cancer to sipuleucel-T. Our mechanistic model incorporates the immune response to the cancer elicited by vaccination, and the effect of androgen depletion therapy. Because only a fraction of patients benefit from sipuleucel-T treatment, inter-individual heterogeneity is clearly crucial. Therefore, we introduce our novel approach, Standing Variations Modeling, which exploits inestimability of model parameters to capture heterogeneity in a deterministic model. We use data from mouse xenograft experiments to infer distributions on parameters critical to tumor growth and to the resultant immune response. Sampling model parameters from these distributions allows us to represent heterogeneity, both at the level of the tumor cells and the individual (mouse) being treated. Our model simulations explain the limited success of sipuleucel-T observed in practice, and predict an optimal combination regime that maximizes predicted efficacy. This approach will generalize to a range of emerging cancer immunotherapies.


Introduction
This year, 1 in 9 men will be diagnosed with prostate cancer (PCa), making this disease the most common cancer diagnosis among American men [1,2]. Despite most cancers seeing a drop in deaths, PCa deaths are on the rise in the U.S. In fact, in 2020, the number of men who died from PCa hit a record, two decade high, with an increase of 5 percent over the previous year. [1]. The risk of death from PCa is highest among men who develop metastatic and/or castrate-resistant disease. PCa also contributes to racial health disparities, as lethality in men of recent African ancestry is roughly twice as high as in men of European ancestry [3]. These statistics are staggering and unquestionably demonstrate the urgency of explaining the limited success of currently used drugs and predicting optimal dose scheduling strategies that maximize therapeutic efficacy.
Typically, the first-line treatment for patients with biochemically failing and metastatic PCa involves blocking the bioavailability of androgens to cancer cells by the constant or periodic application of a combination of chemical castration agents. This is known as androgen deprivation therapy (ADT). Despite initial response rates of 80 to 90 percent, a majority of men treated with ADT progress to a castration resistant state, necessitating alternative interventions [4,5]. Therapeutic advances made over the last decade have conferred varying degrees of survival benefit to patients with metastatic castration-resistant prostate cancer (mCRPC) [6,7]; however, the mortality rate from PCa has first plateaued, and now increased in recent years [1]. Today, mCRPC remains a clinically challenging late-stage cancer with no curative treatment options [7].
The chemotherapeutic agent, docetaxel, was the first therapy to demonstrate a survival benefit for patients with mCRPC [7,8]. In fact, it was the only life-prolonging therapy available from 2004 through 2010 [9]. Recently, the treatment landscape of advanced PCa has evolved and multiple agents including abiraterone, enzalutamide, apalutamide, darolutamide, radium-223, and sipuleucel-T have shown efficacy in improving overall survival, leading to their licensing for the treatment of mCRPC [10,11]. In particular, sipuleucel-T (Provenge) is an autologous cellular immunotherapy approved by the FDA in 2010 for the treatment of patients with mCRPC [6]. Sipuleucel-T is a PCa live cell vaccine consisting of a patient's own peripheral blood mononuclear cells (PBMC). PBMC are harvested via leukapheresis, and cultured in the presence of a recombinant protein, a chimera of prostatic acid phosphatase (PAP) and granulocyte-macrophage colony-stimulating factor (GM-CSF). This chimera activates the mononuclear cells, transforming them into antigen presenting cells (APCs) which are re-injected into the patient 3 days later [6,12].
However, sipuleucel-T has not been widely adopted [13]. This is due, in part, to the limited success it has exhibited clinically. The first phase III clinical trial of sipuleucel-T reported only a modest median overall survival benefit of 4.1 months (25.8 versus 21.7 months over the placebo group), and no significant difference in time to progression (3.7 versus 3.6 months over the placebo group) [14]. Furthermore, the vaccine is given in three doses typically administered every two weeks [6]. However, as the FDA notes [15], an optimal schedule has not been established, with intervals between doses ranging from 1 to 15 weeks in controlled clinical trials.
Our aim is to uncover processes governing the response of cancer cells and individual patients (here, mice) to immunotherapy. From these processes, mathematical models will recapitulate clinical results and predict dosing schedules and therapy combinations with better outcomes. To this end, we develop a nonlinear dynamical systems model of prostate cancer response to sipuleucel-T, alone or in combination with ADT, which incorporates a wide range of biological interactions and processes. We calibrate our model with experimental data on the treatment of PCa xenografts in mice, reported in [16]. Consequently, (virtual) mice are used here as a surrogate for human patients.
We next introduce our novel modeling paradigm, Standing Variations Modeling, which exploits uncertainty in parameter values to derive information on the heterogeneity that is a characteristic feature of cancers. This heterogeneity is key to understanding the results of the sipuleucel-T clinical trial, and informs alternative strategies for maximizing the therapeutic potential of this vaccine. In our approach, we first infer probability distributions from which model parameters most likely arise. This differs from traditional modeling approaches that only use static or mean expression levels and cellular responses, ignoring variation across cell populations or patients. We then sample randomly from these posterior distributions on model parameters, resulting in parameter vectors, each of which can be thought of as a distinct individual within a (virtual) population. The simulated cohort of heterogeneous individuals is then used for treatment optimization and conducting in silico clinical trials. Our Standing Variations Model is inspired by the phenomenon of standing genetic variation in evolutionary biology [17,18]. When selection acts on standing variation, inter-individual (for us, inter-patient or inter-rodent) differences in reproductive fitness arise from existing diversity in the population, rather than from new mutations. In our paradigm, we model divergent trajectories as arising from such existing variation in the population, rather than changes during the trajectory itself. For the models discussed here, this "standing variation", meaning existing diversity in parameter values among the simulated mice, may or may not have a genetic basis.
We build on previously proposed mathematical models of PCa response to ADT [19][20][21][22][23][24][25][26][27][28], as well as models of tumor-immune interactions (we refer the reader to [29,30] for recent reviews). Recent studies have also reported models of PCa response to vaccination alone [31,32] or in combination with ADT [33]. Our approach differs from these in three major respects. Our model of PCa and immune system interactions builds on the tumor-immune interaction models of Radunskaya et al. [34] and Robertson-Tessi et al. [35], and on Rutter and Kuang's model of PCa treatment with ADT and immunotherapy [33] by incorporating a higher level of biological detail, tailored to the specifics of PCa. We will show that these additional relationships drive emergent behaviors important to the outcome of PCa and ADT combination therapy. Moreover, terms in our model are mechanistically derived from biomolecular first principles, whenever possible. Finally, in addition to utilizing experimental data for parameter estimation, our Standing Variations approach reformulates this inverse problem by inferring the distributions from which model parameters arise. Thus, we are able to capture variability within, and across individuals with a deterministic model.

Model Schematic
Our model of PCa growth and response to ADT and immunotherapy is cast as a system of coupled non-linear ordinary differential equations (ODEs) tracking the temporal evolution of key cellular and chemical species. A model schematic is shown in Figure 1.
In our formulation, two physiological compartments are considered where model species reside: the tumor and lymphatic organs such as tumor draining lymph nodes or the spleen. Within the tumor, our model captures androgen-sensitive and castrationresistant cancer cell proliferation, and two forms of cell death, namely, apoptosis and necrosis, which occur in response to environmental stresses or applied treatment [36]. Dead cells are cleared from the tumor site by phagocytes such as macrophages (not included explicitly in our model) and, to a lesser extent, dendritic cells (DCs) [37]. Specifically, dying cells in our model release 'find-me' signaling molecules such as LysoPC and S1P, which trigger an innate immune response resulting in the recruitment of immature DCs to the tumor site [38]. We make the important distinction between cell death via apoptosis and necrosis because immature DCs can only transform into activated antigen presenting cells (APCs) when they phagocytose necrotic tumor cells [39]. Our model also captures APC maturation within the tumor along with their transformation to a phenotype marked by the up-regulation of co-stimulatory molecules such as CD80/CD86, and the expression of various cytokines necessary for the activation of effector T cells [40].
Mature APCs subsequently migrate to lymphoid organs such as tumor draining lymph nodes or the spleen where they activate resident naïve T cells into CD4 + helper T (Th) cells and CD8 + cytotoxic T lymphocytes (CTLs) [40,41]. Activated Th cells, in turn, secrete cytokines such as IL-2 which induces proliferation of all activated T cell populations [42]. Post-activation, our model assumes that all activated T cell populations migrate to the tumor site [41]. To account for fact that the speed of CTL infiltration deceases as tumor volume increases [43], the maximum rate of CTL-induced apoptosis is scaled by the radius of the tumor. Within both, the tumor and the lymphoid compartments, CTLs induce cell death via direct contact in any cell type that expresses their cognate antigen, for instance tumor cells, and activated and mature APCs [44]. We also include regulatory T cells (Tregs) in our framework. Two primary types of Tregs have been identified, namely, naturally occurring or thymus-derived Tregs (tTregs) and peripheral Tregs (pTregs) [45]. tTregs develop in the thymus and may further be classified as resting and activated [46]. The activation of tTregs depends on T cell receptor (TCR) stimulation [45,46], for instance, by mature APCs [47]. On the other hand, pTregs are generated de novo from conventional CD4 + T cells [48], including at the site of the tumor [49]. Indeed, there is evidence that an increase in Treg numbers in cancer patients is a consequence of both, the recruitment of tTregs [49], and the generation of pTregs [48,49]. We account for these processes in our model as follows. For simplicity, we do not distinguish between tTregs and pTregs; rather, resting Tregs are assumed to localize at a constant rate to the lymphoid compartment, where they undergo activation after coming in contact with mature APCs. Additionally, activated Th cells may transform directly into activated Tregs as described below. Once activated, Tregs migrate to the tumor site [50].
Within both, the tumor and the lymphoid compartments, activated Tregs are assumed to exert their immunosuppressive functions in two distinct ways: (i) long range interactions, via the expression of immunosuppressive cytokines such as IL-10 and TGF-β; and (ii) short range interactions involving direct contact with target cells [45,51,52]. Here, we explicitly include TGF-β expression and function, although the model readily generalizes to including additional chemokines. TGF-β induces the conversion of activated CD4 + T cells into activated Treg cells [50,53,54]. This conversion may also occur within the tumor, under signaling from TGF-β secreted by cancer cells, macrophages and activated Tregs [52]. Since we do not model macrophages explicitly, we assume that the source term for macrophage-derived TGF-β is proportional to the numbers of dead cells. TGF-β also exerts an immunosuppressive effect by inhibiting: the production of IL-2 by activated Th cells, the proliferation of activated Th cells and CTLs, and CTL cytotoxicity [55]. Additionally, Tregs down-regulate the activation of the transcription factor NF-κB in DCs, thereby inhibiting their maturation [56]. All of these processes are incorporated in our model. Finally, activated Tregs induce cell death in effector T cells or CTLs in a contact-dependent manner involving granzymes A and B [57].
The effects of ADT and vaccination with sipuleucel-T are incorporated as follows. Under ADT, androgen sensitive cancer cells cease to proliferate and undergo cell death primarily via apoptosis, although some cells may undergo necrosis or necroptosis [58]. Castration resistant cancer cells are assumed to be dependent on androgens to a far lesser extent; consequently, their proliferation rate may slow down under ADT, and limited cell death may be induced. Under vaccination, mature APCs are injected into circulation, from where they may undergo clearance or accumulate in highly vascularized organs such as the spleen [59]. Once in the spleen, the mature APCs are presumed to trigger an anti-tumor immune response as described above.
Further details on model derivation and the full system of ODEs representing the above processes are available in the Supplementary Materials.

Parameter Estimation
Where possible, parameter values were taken from the literature. Values for the remaining parameters were chosen so that model simulations matched available experimental data, in a least squares sense. The results of these fits are shown in Figure 2. In particular, parameters relating to tumor cell proliferation and death, dead cell clearance, immune cell infiltration into the tumor, and chemokine expression were fit to data from PCa xenograft experiments reported in [16]. Briefly, prostate cancer xenografts comprising the Myc-Cap cell line were established in 8-10 week-old male FVB/NJ mice. Tumor diameters were measured periodically, and ADT started on day 22, when tumor volumes reached approximately 400 mm 3 . The onset of castration resistance was observed, and defined as when tumor volumes increased to a minimum of 420 mm 3 after the initial decline post-ADT initiation. This time-course tumor volume data pre-and post-ADT was used to estimate the rates of tumor cell proliferation, the ADT-induced death rate of androgen sensitive cells, and the rate of dead cell clearance. The results of these fits are shown in Figure 2A. Further, Bladou et al. [60] have estimated that roughly 0.66% of cancer cells in prostate cancer xenografts are apoptotic in the absence of any treatment, from which the constitutive or background rates of cancer cell apoptosis were estimated. Soggard et al. [61] have estimated that equal fractions of untreated prostate cancer xenografts are apoptotic and necrotic. Therefore, we take the constitutive or background rates of cancer cell necrosis apoptosis to be equal. The results of these fits are shown in Figure 2B.
In the same set of experiments by Shen et al. [16] described above, the immune presence at the site of the xenograft was also recorded. Specifically, the total numbers of DCs (immature, activated and mature), CD4 + T cells (Th cells and Tregs), CD8 + T Cells, and Tregs were counted prior to ADT initiation, 7 days post ADT initiation, and once castration resistance had emerged. These data were used to estimate the following parameters: immature DC localization rates at the tumor site; activation rates of T cells by mature APCs in the lymphoid compartment; the death rate of activated Tregs; TGF-βmediated rate of Th cell conversion to Treg; and the rates of TGF-β expression by tumor cells and macrophages. The results of these fits are shown in Figure 2C-E.
Further details of this process, including a list of parameter estimates, can be found in the Supplementary Materials.

The Standing Variations Method
Our model has 17 parameters that need to be estimated from the data shown in Figure 2.
Biologically realistic values are assigned to a further 8 parameters, not characterized in the literature. We therefore expect a lack of practical identifiability or estimability in some of these parameters [62]. Furthermore, fixed parameter values, even if subject to uncertainty, would account only for standard error in experimental measurements, and be valid only for an 'averaged' individual mouse. Therefore, the fits shown in Figure 2 reflect samples from a distribution of parameter values which should reflect inter-mouse heterogeneity, and not experimental error. Capturing this heterogeneity is essential in understanding-and predicting-the response of a population to any therapeutic intervention.
Here, we propose the novel paradigm of Standing Variations Modeling, which exploits both, the uncertainty in data, and the inestimability of model parameters, to capture heterogeneity across a population with a deterministic model. Our method is outlined in Figure 3A. We begin by solving the inverse problem of inferring the probability distribu-tions that model parameters follow, rather than estimating their precise values. Specifically, the multivariate uniform distribution is taken as a prior for the unknown model parameters. We then employ sampling importance resampling (SIR), a universally applicable method for obtaining draws from an unknown distribution [63], to infer posterior distributions on these parameters. Draws from the uniform priors are resampled based on their importance ratios, which measure the agreement between the approximated distribution and the experimental data, and are expected to be proportional to the resampling probabilities given the unknown distribution [64]. The SIR method has been successfully applied to quantify parameter uncertainty in epidemiology [65], systems pharmacology [66] and in nonlinear mixed effects models [64]. We selected 18 out of the 25 model parameters where identifiability issues are expected. This choice was motivated based on our biological knowledge of this system, since these parameters are likely to determine PCa response to ADT or immunotherapy. We first used Sobol sequences, a quasi Monte Carlo method [67,68], that allowed us to rapidly converge to a uniform sample from this high dimensional hypercube in parameter space. We then implemented the SIR method, assigning the original samples binary weights, based on whether or not they lie within the standard errors in experimental data points. The complete set of posterior distributions on these 18 parameters is available in the Supplementary Information. For illustration purposes, we plot four representative posterior distributions in Figure 3. As can be seen, the ADT-induced death rate of androgen sensitive cancer cells has narrow almost normal posterior, indicating that this parameter is, in fact, estimable from the given data ( Figure 3B). On the other hand, the posterior on the necrotic fraction of ADT-induced cell death is almost uniform, indicating that this parameter is inestimable from the given data ( Figure 3E).
In the final step of our Standing Variations Modeling approach, these posteriors on model parameters are sampled randomly to generate a collection of parameter vectors, each of whom is representative of a virtual individual. In our case, we generated a cohort of several thousand patients (mice) in silico, each characterized by their individualized parameters. An underlying assumption in this approach is that parameters that are inestimable from the given data, vary widely across the population. Having generated our population of test mice, we can conduct in silico clinical trials designed to predict the response of the population to sipuleucel-T, alone or in combination with ADT, as well as optimize the relative scheduling of these therapies.

In Silico Preclinical Trial
For each preclinical trial simulated here, we first simulated a large enough cohort of virtual mice (as described above) so that we could randomize 5000 mice in each treatment arm. The duration of these preclinical trials was fixed at 100 days. Following [16], the virtual animal was 'euthanized' or presumed dead when its tumor reached a size of 1500 mm 3 . We also assumed the animal was cured if its tumor volume shrank below 1 mm 3 , which is below typical detection thresholds. Overall survival was defined as the time period between treatment initiation and death.

Approximation of Sensitivity Analysis
To assess the contribution of model parameters to the distinct outcomes of a cure (tumor volume < 1 mm 3 ) versus death (tumor volume exceed the survival threshold of 1500 mm 3 ), we utilized a statistical approach: the two categories were fit as an outcome to a multiple logistic regression of the parameter values using the generalized linear model (in R, stats::glm) function in the statistical programming language and environment R. Because the relationship between probability of cure and parameter values may be non-linear, individually significant parameters were then fit to a second multiple logistic regression model as orthogonal polynomials (in R, stats::poly) of up to the third degree.

Treatment Optimization
All treatment optimization was carried out on a cohort of 50,000 simulated mice, using a genetic algorithm [69]. Genetic algorithms mimic the principles of genetic evolution to solve optimization problems. For any desired objective function, such as maximizing the number of survivors or the median survival time, an initial pool of treatment strategies, for instance timing between doses, was created randomly. These strategies were then ranked based on the value of the objective function at the end of the treatment period. The top 20% of strategies were paired to create 'offspring' strategies. These inherit different properties from each of their parents and are also allowed a small chance of mutation. The resultant next generation of treatment protocols were sorted once again based on fitness, and the process repeated until convergence.

Sipuleucel-T Alone Does Not Significantly Improve Overall Survival of Mice
We first investigated the effect on mouse survival, of the two treatments-ADT and vaccination with sipuleucel-T-administered as monotherapies. A cohort of 15,000 virtual mice was generated as described in Section 2.3. A preclinical trial was simulated as described in Section 2.4 with three treatment arms: control, ADT alone or vaccination alone. Our treatment protocols were consistent with those in [16]. That is, PCa xenografts were initiated in each mouse at day 0, and at day 22, when the xenografts reached an approximate volume of 400 mm 3 , the mice were randomly allocated to the treatment groups indicated. ADT was administered continuously to those in the ADT group and, similar to the clinical protocol for administering sipuleucel-T in humans, the animals in the vaccine group were inoculated weekly for a total of three doses.
The resultant Kaplan-Meier survival curves are plotted in Figure 4A. Median survival times for the mice in the control, vaccination, and ADT groups were 30 days, 32 days, and 58 days, respectively, (p-values < 0.0001). The hazard ratio for death in the control group was 4.92 (95% CI, 4.67-5.19) compared to the ADT group, whilst it was 1.18 (95% CI, 1.14-1.23) compared to the vaccination group. Thus, our virtual mice population has a remarkably similar response to vaccination with sipuleucel-T, as the human population, with the vaccine conferring a very modest survival advantage over the control group. Our model suggests several potential reasons for this disappointing result. For instance, the inherent variability in the population, especially in terms of response to the vaccine, coupled with treatment initiation when the xenografts are already quite large (∼400 mm 3 at day 22 post-implantation), and a strong immunosuppressive response from Tregs all contribute to the observations in Figure 4A.
However, a closer inspection of the survival curve for mice on the vaccination arm (blue curve) reveals that seven mice survived to the end of the preclinical trial. Although seven mice out of a cohort of 5000 is small, this finding was a consistent prediction of the model across different cohorts of mice (data not shown). We therefore postulate that in rare cases, individuals may benefit significantly from sipuleucel-T.

Optimization of Vaccine Scheduling
We next investigated whether altering the schedule of vaccination could improve on the median overall survival and/or the number of surviving mice at the end of the trial. The results of this analysis are shown in Figure 4B,C.
Our model suggests that the survival benefit compared to control that vaccination confers is modest, in part, due to initiating treatment when the tumor is already large. We therefore varied the time of the first dose of vaccine on a fixed cohort of 5000 mice. The earliest a dose could be given was at day 11 post-implantation, when the xenografts became palpable [16]. Clinically, this mimics starting vaccination at the time of first diagnosis. The remaining two doses were given weekly, as before. Figure 4B plots a graph of the number of surviving mice at the end of this trial, as a function of days post-implantation when vaccination was initiated. Starting treatment as early as possible confers a significant survival benefit to the population of mice, with the number of survivors increasing to 50 when vaccination was initiated on day 11 as compared to just 7 when vaccination was administered on day 22 ( Figure 4B, dark shaded square). This result is explained by the fact that in our model, we account for a slow-down in the speed of CTL infiltration into the tumor with an increase in tumor volume [43]. Figure 4C shows the difference in survival/failure times for control mice, and for mice in vaccination schedules. The mice are divided into quartiles independently, to highlight heterogeneous benefits of early vaccination. Roughly 25% of mice (the first quartile) see no benefit from early vaccination. Another 25%, the fourth quartile, corresponding to mice with a relatively long survival even under control, show a 5 day increase in survival time for early vaccination. None of the simulated mice benefit appreciably from late vaccination, at 22 days.
An optimal dosing schedule has not been established for treating human patients with sipuleucel-T [15]. We therefore sought to maximize the median survival time by varying the intervals between successive doses of the vaccine. A cohort of 50,000 virtual mice was generated as described in Section 2.3. The timing of all three doses of vaccine were allowed to vary independently, with inter-dose times varying between 1 and 21 days. The earliest that treatment could be applied was at day 11 post-implantation (see comment in preceding paragraph). The optimal dosing schedule was identified by employing a genetic algorithm as described in Section 2.6. Unsurprisingly, initiating treatment at the earliest possible timepoint (day 11) was optimal. Figure 4D plots the median survival time (x-axis) as a function of inter-dose intervals. The color of the dots indicates the interval between the first and second vaccine dose with blue dots corresponding to longer waiting times, while the size of the dots indicates the interval between the second and third vaccine dose with larger dots corresponding to longer waiting times. Median survival time was maximized when all three doses were administered 3-4 days apart. However, the median survival time only improved by 1.4 days, going from 32 to 33.4 days on this optimal protocol. These results offer an explanation as to why an optimal dosing schedule for sipuleucel-T has not been established in the clinic. Varying the inter-dose intervals across a large range of values only yielded a modest improvement in survival, which is possibly what has been observed in clinical trials of sipuleucel-T.
Encouraged by our finding that early administration of vaccine improves the survival chances of mice, we instead sought to maximize the number of survivors at the end of the trial, on the same cohort of 50,000 mice. Once again, initiating treatment at the earliest possible time-point (day 11) was found to be optimal. The number of surviving mice are plotted on the y-axis as a function of inter-dose intervals, in Figure 4D. A maximum of 580 mice (1.16%) survive, of which 395 are cured (tumor size < 1 mm 3 ), if the first two vaccine doses are given within 3 days of each other, whilst the third dose is administered 13-15 days later. This is a marked improvement from the 7 survivors out of 5000 (0.14%) when mice were vaccinated starting on day 22, on a weekly schedule. We remark that the number of surviving mice (equivalently, the higher quantiles of the survival time) is more sensitive to differences in strategy than is the median survival time. This is typical in biomedical studies where only a subpopulation of patients benefit from the intervention (for a review, see [70]).

Characteristics of Response to Vaccine
The characteristics of the surviving mice-and in particular, the rare, cured subpopulation-may be of crucial importance in determining who will benefit the most from vaccination, and to what extent. The power of a mechanistic model is to extrapolate into these desirable regions of parameter space which can then be targeted and explored in a clinical setting.
We therefore sought to identify those parameters that are significantly associated with a cure (tumor size < 1 mm 3 ). Table 1 gives trends in model parameters significantly associated with this outcome in an individual simulated mouse, from the final cohort of 50,000 mice. Multiple logistic regression compared 395 cured mice to 49,502 'euthanized' or dead mice, with mice surviving-but-not-cured at the end of the simulation excluded from this analysis. Positive effect sizes indicate terms of polynomials of parameters that are higher in the cured mice, while negative effect sizes indicate the opposite. The p-values can be made arbitrarily small by increasing the number of simulated mice, but the relative value of Z-scores between parameters reflects their contribution to the outcome, approximating a sensitivity analysis if the model could be analyzed in such a way. Non-significant parameters are not reported in the table, for the sake of brevity. Table 1. For each parameter, including higher-order terms where significant: the Dimensions of the parameter, the Effect Size (on a logistic scale, per unit), the standard error in the estimate of the effect size (SE), the corresponding Z score and p-value from the multiple logistic regression.

Dimensions
Effect The significant parameters were found to be: v tum (10), the xenograft tumor volume at day 10 (1 day prior to vaccination initiation); α N , the rate of cancer cell proliferation (since we are administering vaccination alone, we do not distinguish between an androgensensitive versus castration-resistant phenotype); δ D , the rate of dead cell clearance by (assumed) macrophages; λ T8 , the maximum rate of naïve CTL activation by mature APCs; α C T , the rate of TGF-β expression by cancer cells; α D T , the rate of TGF-β expression by (assumed) macrophages, taken here to be proportional to the number of dead cells; δ Tr , the rate of activated Treg death; K T , the half-saturation constant for APC-mediated activation of T cells; and K I , the CTL to target cell ratio at which rate of CTL-induced cell kill is half its maximum value. For most parameters, first and second order terms have opposite effects, indicating "diminishing returns"; exceptions being α C T , which shows a linear relationship, and K I , which shows a compounding effect. The overall differences in parameter values are summarized as box plots in Figure 5. Parameters are normalized to mean and standard deviation across all simulated mice. Only those parameters are represented that were found to be significantly associated with a cure (tumor size 1 mm 3 ) when vaccination is administered as a monotherapy, on the optimal schedule predicted by the model. Red boxes correspond to dead mice, and teal boxes correspond to cured mice.
From the Z-scores in Table 1, we conclude that the most important determinants of response to vaccination are α N and δ Tr , with a low tumor cell proliferation rate or a high rate of activated Treg death associated with cured mice. The next largest effects come from δ D , α C T , α D T and K T , with high values of δ D and low values of α C T , α D T and K T associated with cured mice. δ D , α C T and α D T characterize the amount of TGF-β in the system: a large δ D would result in more efficient clearance of dead cells from the tumor space, and thus correspond to reduced macrophage numbers, and a smaller source of TGF-β. Equally, small α C T or α D T would indicate diminished TGF-β production and consequently, elevated CTL activity.

ADT and Vaccination Combination Therapy
We investigated the potential of combining sipuleucel-T with ADT, and sought an optimal relative dosing schedule that maximizes median overall survival. As before, optimization was carried out on a cohort of 50,000 simulated mice, with timing of ADT initiation fixed at day 21 post-implantation. Three doses of vaccine could be administered with inter-dose timing allowed to vary between 1 and 21 days. The earliest that vaccination could be started was at day 11 post-implantation. The optimal dosing schedule, identified by employing a genetic algorithm, was vaccine doses at days 12, 16 and 27 post-implantation. That is, our model suggests administering two doses of vaccine in quick succession, and waiting for the third dose until after ADT initiation.
To illustrate the potential benefit of combining vaccination with ADT across a population, we conducted an in silico preclinical trial, comparing: (1) our predicted optimal dose schedule; (2) initiating vaccination simultaneously with ADT (vaccine given in three weekly doses); and (3) administering vaccination 7 days post-ADT, when castration resistance is observed in the mice xenografts [16]. We remark that this last schedule mimics the current clinical protocol for administering sipuleucel-T. As before, 15,000 virtual mice enrolled in our trial and 5000 mice randomized into each treatment arm. The resultant Kaplan-Meier survival curves are plotted in Figure 6A. Median survival times for the mice in the clinical arm (schedule (3)), simultaneous arm (schedule (2)), and optimal arm (schedule (1)) were 58 days, 62 days, and 64 days, respectively, (p-values < 0.0001). The hazard ratio in the clinical group was 1.58 (95% CI, 1.52-1.65) compared to the simultaneous group, whilst it was 1.79 (95% CI, 1.72-1.87) compared to the optimal schedule group. Moreover, only 20 mice (0.4%) were predicted to survive to the end of the trial on the clinical protocol. When sipuleucel-T was co-administered with ADT, there were 284 survivors (∼5.7%) with 179 cured, whilst this number went up to 441 (∼8.8%) with 318 cured on the optimal protocol. Thus, our model suggests that instead of waiting until after ADT has failed, vaccination should be given prior to starting androgen deprivation, and in a staggered, rather than periodic fashion. This is predicted to reduce the rate of cancer death by approximately 45%.
Finally, we sought to identify those parameters that are significantly associated with a cure (tumor size < 1 mm 3 ), when ADT is combined with vaccination and administered per the optimal schedule found above. In addition to the parameters found to be significantly associated with a cure when vaccination is administered as a monotherapy (see Section 3.3), now α M , the rate of castration-resistant cell proliferation, and M(10), the number of castration resistant cells that are present in the tumor prior to treatment initiation are also predicted to be significant. The overall differences in parameter values are summarized as box plots in Figure 6B (see also  Only those parameters are represented that were found to be significantly associated with a cure (tumor size 1 mm 3 ) when ADT and vaccination are administered as a combination, on the optimal schedule predicted by the model. Red boxes correspond to dead mice, and teal boxes correspond to cured mice.

Discussion
The landscape for PCa therapy has developed rapidly with the approval of new classes of treatments, including immune-based approaches, for use in various PCa disease states [6]. FDA approval of sipuleucel-T for the treatment of castration resistant PCa marked the first cancer "vaccine" approved for use in a treatment setting and sparked renewed interest in cancer vaccines more broadly [71]. Despite the initial enthusiasm for sipuleucel-T, it has not been widely used compared with other available therapies, in part, due to its limited clinical success [14,72].
A powerful and practical way to analyze potential mechanisms driving therapeutic success and failure, and to optimize dosing schedules of novel drugs, is to use multi-scale computational (dynamical systems) modeling of cancer response to treatment. Biomedically sound parametrization is essential to the validity of such modeling approaches. Validated mechanistic models enable reliable, quantitative predictions [62,73]. However, data are typically sparse or incomplete, leading to parameter identifiability issues [62]. Furthermore, even if some parameter values can be determined precisely, these would only be valid for an 'averaged' individual (cell/mouse/human), and cannot capture the inherent variability across a population, let alone within an individual.
In this paper, we developed a detailed and integrative mathematical model of PCa response to sipuleucel-T and ADT, in order to shed light on the limited clinical success of the sipuleucel-T vaccine. Our model simulates the growth and response to treatment of PCa xenografts in an individual (mouse). However, we aim to understand the response of a population to vaccine therapy. Therefore, we proposed a novel modeling paradigm that exploits uncertainty and variability in data to capture heterogeneity across a population. This heterogeneity is key to understanding both health disparities in PCa outcomes (including especially ethnic/racial disparities [3]), and the results sipuleucel-T clinical trials. We term this synthesis of importance sampling and deterministic differential equation modeling a 'Standing Variations Model', inspired by similar approaches in population genetics [17,18].
Using our modeling framework, we first conducted an in silico preclinical trial of sipuleucel-T administered as a monotherapy. Survival data from our simulated cohort of mice was found to be in good qualitative agreement with results from the first phase III clinical trial of sipuleucel-T [14], with vaccination providing only a modest survival benefit compared to control. In our model, several parameters that are key determinants of response to vaccination vary significantly across the population, and no one parameter in and of itself determines a favorable outcome. However, these parameter differences do suggest mechanistic causes of limited sipuleucel-T efficacy, in both our simulated mice and in the clinic. For example, tumor sizes are likely too large at the time of treatment initiation, making it almost impossible for active CTLs to successfully infiltrate the tumor and eliminate target cancer cells. Further, the vaccine elicits a strong immunosuppressive response in terms of Treg activation and TGF-β expression within the tumor, both of which appear to limit efficacy.
We next varied vaccination schedules to: (1) understand why an optimal dosing schedule has not been established for sipuleucel-T [15]; and (2) what such a schedule would be. Our simulations indicated that even varying the time interval between vaccine doses over a wide range could not appreciably improve on median overall survival of the mouse population. This is again in agreement with what has been observed clinically for sipuleucel-T [15]. However, a small fraction of the population was predicted to benefit significantly from vaccination, motivating us to instead seek treatment protocols that maximize the number of such individuals. A non-standard dosing schedule-two quick doses followed by a longer interval to the third dose-was predicted to be optimal. These results support an off-label trial, with investigational use of sipuleucel-T in patients with androgen-sensitive cancer to rigorously assess potential survival benefits.
Sipuleucel-T is currently approved for the treatment of patients with mCRPC, that is, after ADT has failed [6]. We therefore conducted a preclinical trial to investigate the potential of combining ADT with vaccination. Our model predicted administering 2 doses of the vaccine prior to ADT was optimal, in terms of maximizing the median overall survival, and resulted in 6.36% cured mice as compared to 0.08% on the schedule that mimics the current clinical protocol. This suggests that there is potential to significantly improve the clinical performance of sipuleucel-T by combining it with other therapeutic modalities such as ADT. Furthermore, a multiple logistic regression analysis was conducted to identify the parameters that are significantly associated with a cure (tumor volume < 1 mm 3 ) under vaccination alone, or vaccination in combination with ADT. Parameters associated with a substantial immunosuppressive response to the vaccine, such as TGF-β production rates, were strongly negatively correlated with this favorable outcome. These results highlight a potential for combining anti-TGF-β therapy with vaccination in improving outcomes, although this would require investigational agents such as Galunisertib [74].
The standing variation modeling approach introduced here is flexible enough to include these agents, as well as to model adaptive clinical trial designs incorporating any of these parameters as biomarkers. Work is ongoing on future iterations of this model which will include additional combination therapies and design elements, as well as additional mechanistic features such as the activation or differentiation of monocyte lineage cells (for instance, macrophages) and incorporating immunosuppressive functions of cancer cells via the PD1-PDL1 axis [75].
Clinical translation of biological insights, such as those underlying cancer vaccines, requires tremendous investments of time, money, and limited clinical-scientific resources such as trained physician-scientists and the recruitment of rare patient populations. Mathematical modeling can both, predict the success of these investments, and optimize that success given that even safe and effective novel therapeutic modalities may provide limited clinical benefit if administered in a sub-optimal way. This optimization is a particular strength of our approach, since clinical trials randomizing hundreds of thousands of mice or patients to different combination therapies, each administered over a range of potential schedules, would be prohibitively expensive.  Figure S2.1: Posterior distributions of parameters relating to immune presence in the tumor inferred using a sample importance resampling algorithm. Figure S2.2: Posterior distributions of parameters relating to tumor growth and response to ADT inferred using a sample importance resampling algorithm. Figure S2.3: Assumed distribution of K I , the CTL to target cell ratio at which rate of CTL-induced cell kill is half its maximum value. Funding: This work was partially supported by NIH/NCI U01CA243075 (TLJ).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments:
The authors thank Krishna Natarajan and Giray Okten for many helpful discussions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: