Myeloid-Derived Suppressor–Cell Dynamics Control Outcomes in the Metastatic Niche

Myeloid-derived suppressor cells (MDSC) are immature myeloid immune cells that become pathologically activated with potent immunosuppressive activity (1–7). Since the introduction of the term “MDSC” in the late 1990s (4–6), there has been a great deal of effort to understand MDSC phenotypes and dynamics. MDSCs are implicated in the regulation of immune responses in many biological contexts and pathologic conditions, including cancer, inflammation, wound healing, and autoimmune disorders (1). Some have gone as far as to claim that MDSCs are “the most important cell you have never heard of” (8). Recently, with the advent of high-dimensional measurement technologies, including mass cytometry and single-cell RNA sequencing, the characterization of MDSCs and their roles in diverse contexts has become more refined (7, 9). Here, we characterize MDSCs by their function—immunosuppressive activity—rather than their expression phenotype (e.g., CD11b⁺ and Gr-1⁺ in mice), bypassing the need to delve into the heterogeneity of the CD11b⁺Gr-1⁺ population at the single-cell level.

In the context of cancer, the role of MDSCs is convoluted, in part due to the complexity of the tumor microenvironment and related immunology (3, 10–14). MDSCs certainly play important roles in tumor microenvironments (8, 9, 15, 16); increased levels of MDSCs are associated with poor clinical outcomes (2, 12, 17–19); although an important caveat is that studies often measure only circulating MDSCs. There is compelling evidence that MDSCs can effectively shield tumors from antitumor immune responses mediated by cytotoxic T cells and natural killer (NK) cells (20–24). Thus, targeting MDSCs as a way to sensitize nonimmunogenic tumors is an attractive treatment strategy in cancer immunotherapy (16, 17). MDSC dynamics have also been studied in the specific context of breast cancer, where they have been shown to affect the progression of primary breast tumors and associated metastases (7, 15, 18, 23, 25–27).

Understanding tumor–immune–MDSC dynamics is by nature a systems biology problem. Mathematical and computational modeling are essential to tease apart the intricate relationships involved (28, 29). There have been relatively few works, in comparison to experimental/clinical interest, in the literature that report mathematical models of MDSCs (30–33). Shariatpanahi and colleagues (30) developed a model described by ordinary differential equations with which they explored therapeutic strategies that aim to restore antitumor immunity, in comparison with experimental data (23). Allahverdy and colleagues (31) developed a stochastic agent-based model that was used to explore the effects of different drugs on MDSC and tumor dynamics. Liao and colleagues (32, 33) developed a model described by partial differential equations that was used to determine optimized drug treatment and to understand primary drug resistance. While these models offer insight into the roles of MDSCs, a rigorous treatment of MDSC dynamics in the tumor microenvironment, fitting models to data, and taking into account the effects of noise remains lacking.

Here, we focus on the effect of MDSC dynamics on metastatic tumor growth following an initial seeding event. Most cancer deaths are a result of metastasis (34): a highly dynamic and stochastic process. Most metastatic tumors are seeded by a small number of circulating tumor cells (13, 34). MDSC migration to the site of a new tumor has been identified as crucial for cancer progression, both in primary tumors and metastases, but the interactions involved are not well understood, in part due to a relative lack of MDSC characterization, the complex tumor–immune environment, and the difficulties associated with tracking cell–cell interactions in vivo (13, 35, 36). As a result, there are many open and pressing questions regarding MDSCs and tumor metastasis (37). How much therapeutic benefit can be gained by blocking MDSC recruitment to the tumor site? Would therapies that decrease the number of circulating MDSCs achieve similar or greater effects? There are now various methods to target MDSCs in peripheral lymphoid organs and their migration to tumor sites. However, it is not clear whether either of these methods alone will be sufficient to inhibit MDSC immunosuppressiveness at a tumor site or whether combination approaches will be required.

To address these questions, we developed a stochastic delay differential equation model of metastatic tumor growth. We included an MDSC delay that could represent delays in MDSC recruitment to the metastatic tumor site as well as delays in MDSC activation to suppress antitumor immune cells. Stochasticity was included due to the inherent noise in the cell dynamics, and to allow the assessment of the probabilistic events of new metastases. We first demonstrated the importance of MDSCs in the tumor immune microenvironment and established conditions necessary for metastatic growth for the deterministic model. We then identified the most important parameters and interactions in the system, to shed light on the underlying biological dynamics. Next, through simulation we explored the impact of MDSC delays on metastatic growth, finding that under certain conditions, inhibiting MDSC recruitment alone might be a highly effective treatment strategy. Finally, we performed Bayesian parameter estimation of models fit to individual tumors growing in vivo, from which we determined tumor- and MDSC-specific parameters. Inference results revealed that knowledge of MDSC-specific parameters is important to be able to accurately predict metastatic outcomes.

Mathematical modeling of tumor–immune cell interactions has been increasingly recognized as critical for understanding strategies to mount an effective response to cancer initiation, spread, and evolution (28, 29, 38–43). Herein, we first describe a theoretical basis for MDSC dynamics in the context of a metastasizing tumor (e.g., in the lung, bone, or liver; ref. 44) from a primary tumor in the breast. For parameterization of the model, we focused on the lung, as it is one of the most common distant metastases sites of breast cancer (45). Thus, we primarily considered an NK cell–abundant environment (46).

Our mathematical model is comprised of four nonspatial delay differential equations to describe tumor–immune interactions incorporating MDSCs. Our model is built on previous literature in mathematical oncology (30, 40) with a focus on the role of MDSCs in metastases (37). We used our previous preclinical data (15), and data reported in the literature (cited below and in Table 1), and focused on important interactions between tumor, immune and MDSC populations, leading to a relatively simple model that allowed us to gain insight into system dynamics and metastatic tumor spread. We included the antitumor immune populations of cytotoxic T cells (CTL) and NK cells. MDSC–CTL interactions are important given the primary function of intratumoral MDSCs is suppression of CTLs (1, 6, 15, 16). MDSC–NK cell interactions are also important (20–22, 24, 25), and NK cells are increasingly being studied as an immune population specifically affected by tumor cells to promote metastasis (47, 48). A schematic diagram of the model is provided in Fig. 1A.

We denote x_T, x_MDSC, x_NK, and x_CTL as the populations of tumor cells, MDSCs, NK cells, and CTL cells, respectively, at time t. The model derived can be expressed conceptually (i.e., agnostic as yet to the form of the dynamics) as follows, where δx_i denotes the rate of change of x_i, i ∈ [T, MDSC, NK, CTL].

where ξ_i(t) is the size of the i^th population, that is, we assumed multiplicative noise (49, 50). We studied the tumor–immune dynamics under the assumption of multiplicative noise given the mounting evidence that biological systems more often exhibit dynamics generated from multiplicative noise models (51).

For the deterministic dynamics, we have:

with description of the parameters in Table 1. We modeled tumor growth [first term of Eq. (C)] according to a Gompertzian model (30, 38), with maximum size η. Tumor cells can be eradicated by NK cells and CTLs (antitumor response), with rates β₁ and β₂, respectively. MDSCs are activated due to their basal circulation, α₂, and die at rate ζ₂. In the presence of tumor cells, immune-suppressive signals (α₃) lead to increased MDSC production, activation, and recruitment to the site of the tumor, and the size of the MDSC population increases; that is, tumor growth can dramatically increase the number of MDSCs. In the absence of tumor cells, α₃ is zero and the MDSC population is unaffected.

MDSCs are produced primarily in the bone marrow, from which they migrate to peripheral lymphoid organs and tumor sites in tumor bearing hosts (13, 52). The delay in activation/recruitment of MDSCs to the site of a tumor was modeled using a Mackey–Glass delay term (53), with a delay of τ [second term of Eq. (C): second row]. We also considered an alternative form for the delay, for which we found the dynamics were similar (see Supplementary Text Section S1; Supplementary Fig. S1). We considered delays only in x_MDSC; while other delays, for example, in CTL activation, may well be important in some contexts, here such a delay was observed to exert only minor effects on the tumor dynamics, due to the low circulating levels of CTL cells (Supplementary Text Section S2; Supplementary Fig. S2). We note that the model does not include MDSC subtypes or maturation, but accounts only for their functional significance as immature myeloid cells with immunosuppressive capability. Future work could include MDSC maturation into other cell types as influenced by the tumor microenvironment (see the Discussion for further details).

Regarding the antitumor immune dynamics, NK cells are produced at rate α₄, and CTLs are activated by the NK cell–tumor cell interaction at rate α₆. In line with (30), both NK cells and CTLs can be activated by the tumor (at rates α₅ and α₇, respectively). We assumed that NK cells and CTLs can be inhibited by MDSCs (at rates β₃ and β₄, respectively), and are lost due to cell death (at rates ζ₃ and ζ₄, respectively).

In simulations of new metastases Eq. (C), the initial conditions are set by the tumor-free steady state Eq. (F), except that we seeded tumor growth by one or two initial tumor cells. Unless explicitly stated otherwise, all parameter values used for simulation were as defined in Table 1. Standard error is defined as standard deviation/|$\sqrt {{\rm{number\ of\ simulations}}}$|⁠.

In the work below we considered analyses of the full SDDE model, as well as various reduced models. In the case that g = 0, the SDDE model reduces to a deterministic delay differential equation (DDE) model. In the case that g = 0 and τ = 0, the model reduces to an ordinary differential equation (ODE) model. All models were developed in the Julia programming language (54), using DifferentialEquations.jl (55). For simulation of the full model, we use the SOSRI algorithm for stiff stochastic differential equations (56). Metaprogramming in Julia enables transitioning between model formulations (SDDE, DDE, or ODE) with ease (57).

Parameter sensitivity analysis was performed to assess the relative importance of parameters on the model given by Eq. (C). We used Morris global sensitivity analysis (GSA; refs. 58, 59), with the steady state of the tumor population defined as the output variable. Table 1 contains the ranges used for GSA as well as the parameter descriptions. Hyperparameters used for the Morris algorithm (implemented in DifferentialEquations.jl; ref. 55) were total_num_trajectory = 1000 and num_trajectory = 100.

RECIST criteria have been developed for use in clinical trials to determine the change in tumor burden of selected target lesions to inform whether a patient is responding to a given therapy (60). We implemented Bayesian parameter inference to fit the model to tumor responses using RECIST to classify tumor sizes and responses over time (61). We fit differential equation-based models to RECIST data following a similar conceptual framework to (38). In the case of our model, we also fit certain MDSC parameters, such as the interaction strengths between the MDSCs and other immune/tumor populations, to assess the effect of MDSC dynamics on clinically-relevant tumor growth. We employed Bayesian parameter inference (62) implemented in Turing.jl (63).

We used in vivo tumor data from a study evaluating the efficacy and safety of anti–programmed death ligand 1 (PD-L1) atezolizumab in advanced non–small cell lung cancer (61). These data were also recently used to fit mathematical models of tumor growth (Study 1, ref. 38). Each tumor was assessed at baseline, before the initiation of treatment in the clinical trial (for the purposes of fitting we set the time of the baseline assessment to be zero). Tumor size was then reassessed approximately every six weeks for 12 months, then every nine weeks, and then at disease progression. At each assessment the tumor size was measured in millimeters in one dimension (x), which we converted to a volume following the convention adopted by Laleh and colleagues, that is, taking the volume (mm³) as |$\frac{1}{2}{x}^3\$|(38, 64). We estimated the number of tumor cells from this volume by multiplying by a factor of 10⁷ (65). From the available data we selected six measurable tumors from six different patients that each have data from at least five time points (including the baseline assessment), were from all three study cohorts, and were representative of the range of the dataset (i.e., tumors that increase/decrease at a variety of rates). We fit the relative change in the tumor population, which is measured as the difference between the measurement and the baseline assessment, divided by the baseline assessment [(measurement-baseline)/baseline, which produces a real number ∈ [−1,∞)]. As the relative change at the baseline assessment was always zero, we removed this data point for all tumors. Since only the tumor data were available, we fit the log transformed data from this population (i.e., log(x_T +1)). All of the data for each of the six tumors is available in the supplementary file tumor_data.xlsx.

For inference, a three-dimensional parameter space was considered, in which we fit the following parameters: β₃ (NK cells inhibition rate by MDSCs), α₆ (CTL stimulation by tumor–NK cell interaction), and α₁ (tumor growth rate). We did not infer all parameters simultaneously due to the high dimensionality of the parameter space; we chose parameters to infer based on their importance as judged by prior knowledge and global sensitivity analysis. The parameters chosen affect all populations of the model. As no information on time since incidence was available, we set the initial conditions according to previous simulations (see Fig. 1B and Table 1) at day 100 (x_T(0) = 8,395.4, x_MDSC(0) = 804.1, x_NK(0) = 197,565.7 and x_CTL(0) = 1,654.4). Therefore, we rescaled η = 10⁵ (tumor maximum size), and all other parameters were set to be as in Table 1 with τ = 0.

A prior distribution over the three-dimensional parameter space was selected using weakly informative normally distributed priors, with means set to the values given in Table 1, and relatively large variances:

where σ is the observational noise. For each tumor, we ran four independent Markov chain Monte Carlo (MCMC) simulations with 2 × 10³ iterations using the No U-Turn Sampler (NUTS) with a target acceptance ratio of 0.65 (66).

In addition, we fit the model to the same six tumors in a CTL-abundant environment. Here we assumed the initial conditions were x_T(0) = 16,654.8, x_MDSC(0) = 707.8, x_NK(0) = 16,162.1 and x_CTL(0) = 778,128.8 (simulation values at day ten) and that the prior distributions for the parameters were:

We used a decision tree classifier to test whether different combinations of the marginal posterior parameters obtained from Bayesian parameter inference could classify tumor responses as either decreasing or increasing over time. Decision trees were built using DecisionTree.jl (67) and cross validation was performed using scikit-learn (68).

As described above, the primary tumor microenvironment of focus is that of the lung, an NK cell–rich environment (46). However, we also considered tumor–MDSC dynamics in a microenvironment that—relative to the lung—is low in NK cells and abundant of CTLs. We modeled this CTL-abundant environment by adjusting key CTL/NK cell–related model parameters (see Table 1 for parameter descriptions). To model a CTL-abundant environment, we considered the parameter values given in Table 1, but decreased the NK cells’ circulating rate to α₄ = 10³, and increased the CTL stimulation by tumor/NK cell–interaction rate to α₆ = 2 × 10⁻³. This led to an increase in the number of CTL cells and a concordant decrease in the number of NK cells.

All code and data are available at a public github repository located here: https://github.com/maclean-lab/ModelingMDSCs. Tumor data (38, 61) are available in the supplementary file tumor_data.xlsx. All other data are available in the figures, tables, and Supplementary materials or from the corresponding author upon reasonable request.

We studied MDSC dynamics in the context of a metastasizing tumor, specifically we focused on breast-to-lung metastasis, that is, metastatic growth in the lung resulting from a tumor in the breast. To parameterize the model (model description in Materials and Methods), we considered the immune cell composition known to be present in tumors in the lung, which is an environment abundant in NK cells (ref. 46; Fig. 1A). We began by analyzing the behavior of the deterministic model [delay differential equations (DDEs); Eq. (C), g = 0]. Simulation of the DDE model for different sizes of MDSC delay (τ) showed that the delay in the recruitment of MDSCs to the tumor site played a critical role in determining metastatic tumor size after one year (Fig. 1B–E). Tumor growth resulted in large increases in the MDSC population. We also saw that increasing τ led to slower growth and smaller sizes of both the MDSC and tumor populations. Increasing the delay led to a lag before the MDSCs received activation signals from the tumor and began to proliferate. Smaller MDSC population sizes led to slower growth/smaller tumor population sizes because a smaller MDSC population made the tumor more susceptible to killing by NK-cell and CTL populations. Note that, given the parameters in Table 1, the same steady state would be reached for any finite τ, 0 ≤ τ < ∞. The time until steady state was positively correlated with the delay τ.

In the case of no tumor (x_T = 0), the tumor-free fixed point of the model is:

If we assume that the tumor reaches its carrying capacity, η, then the tumor steady state is given by Eq. (G) with |$\hat{x}$|_T = η. A comparison of simulations with different carrying capacities is given in Supplementary Text Section S4 and Supplementary Figs S3 and S4.

Using the model, we can also determine whether a small initial number of tumor cells will grow to a nonzero positive steady state (for simulation, we considered steady states with x_T> 10) or will decay. This is a helpful characterization, as metastases are thought to be seeded from a small initial number of circulating tumor cells (1, 13, 34, 37, 69). When we began at the tumor-free steady state Eq. (F), and increased the number of tumor cells by one or two, then taking the highest order terms in Eq. (C) (first line) we saw that the rate of change of the tumor population was initially positive if G > 0. Here G can be defined as the tumor growth threshold, or equivalently, the tumor basic reproductive ratio (analogous to R₀ in epidemiological models; see Supplementary Text Section S5 and Supplementary Fig. S5 for details). G is given by:

Examples of simulations starting from the tumor-free steady state Eq. (F) but with the addition of a single tumor cell are shown in Fig. 2A-C. In these examples, the tumor population grows initially if and only if G > 0. In Fig. 2A the parameter values are as defined in Table 1, giving G ≈ 0.8, and a resulting tumor size at steady state of 9.8×10⁶. We changed the tumor cell death rate (ζ₁) to vary G: to G ≈ 0 (giving a tumor steady state of ≈ 1; Fig. 2B), and to G ≈ −0.2 (giving a tumor steady state of < 1; Fig. 2C). The threshold G thus gives an approximation of whether small numbers of tumor cells will grow into fully developed metastases, of relevance for cancer prognosis, treatment, and progression (36).

We performed parameter sensitivity analysis to assess the relative importance of model parameters on the growth and final size of the tumor population (see Materials and Methods for description of the model and the parameter sensitivity analysis). Because the tumor steady state is independent of the MDSC delay as t → ∞, for sensitivity analysis we set the delay τ = 0.

As seen in the model Eq. (C), the MDSC-specific parameters are α₂, α₃, ζ₂, β₃, β₄, and γ₁. The Morris global sensitivity values of parameters on the (numerical) tumor steady state are given in Fig. 2D, where the MDSC-specific parameters are marked by hexagons. The green (red) color denotes parameters that were positively (negatively) correlated with the tumor steady state. We found that ζ₂ (death rate of MDSCs) and γ₁ (steepness of MDSC production) were the only MDSC-specific parameters negatively correlated with the tumor population, as increasing either of these parameters resulted in fewer MDSCs and thus a more immunosusceptible tumor population.

As shown in Fig. 2D, β₃ (inhibition of NK cells by MDSCs) was the most important MDSC-specific parameter for the tumor steady state. This is in part due to the initially large size of the NK-cell population [ref. 46; see NK cells in Fig. 1B–E and Eq. (F)]: the MDSC population must effectively suppress NK cells for the tumor to be able to grow rather than die out. Similarly, β₄ (inhibition of CTL cells by MDSCs) was important, but less so than β₃ as the CTL population is initially small thus less important for initial growth of the tumor [see CTL cells in Fig. 1B–E and Eq. (F)]. For consideration of the model in a CTL-abundant environment, see the final section of the Results. Figure 2E and F show the effects of parameter β₃ (inhibition of NK cells by MDSCs) on the tumor steady state at both ends of the range studied. We saw that small values of β₃ (Fig. 2E) resulted in small tumors (tumor steady state approx. 2.5×10²), whereas large values of β₃ (Fig. 2F) resulted in much larger tumors (tumor steady state approx. 9.9 × 10⁶).

Comparison of MDSC-specific parameters (hexagons) with non–MDSC-specific parameters (circles, Fig. 2D) showed that α₂, α₃, α₁, η, β₃, ζ₃, β₄, ζ₄, γ₂, and γ₃ were positively correlated with the tumor population steady state; all other parameters were negatively correlated. The most important parameters (as measured by their effect on the tumor steady state) were α₆, β₁, β₂, β₃, and β₄. α₆ is the rate of CTL stimulation by tumor–NK cell interaction, β₁ and β₂ are inhibition rates of tumor cells by NK and CTL cells, and β₃ and β₄ are inhibition rates of NK and CTL cells by MDSCs (see Table 1). Therefore, model dynamics are predominantly influenced by inhibition/stimulation between the competing populations (see Fig. 1A for schematic). This reinforces recent literature that highlights these tumor–immune interactions (especially in the context of MDSCs) as important determinants of outcomes (1, 8, 10, 17, 20, 21, 42).

We next analyzed the stochastic dynamics of the model. Given the seeding of metastases by one or a few cells, stochastic effects are likely to play a large role in the system. To study metastatic tumor establishment and viability we simulated the SDDE model Eq. (C), with MDSC delay τ ≥ 0.

Stochastic simulations allowed for the probabilistic analysis of “successful metastases.” In the deterministic setting, using the parameters defined in Table 1, a metastatic tumor is always formed (G > 0). In the stochastic setting, this is no longer the case. Model outcomes vary even for identical initial conditions due to noise in the system (10, 70, 71). Although we did not study the sources of biological noise here, we expect the major component to result from noise in the intercellular signaling processes, that is, extrinsic noise (72).

To study the probability that a small number of pioneering cells will establish a new metastasis, we implemented simulations that started with two tumor cells (see Supplementary Text Section S6; Supplementary Fig. S6 for larger initial conditions), and denoted a metastasis successful if the number of tumor cells did not drop below one (i.e., |$\lfloor {x}_T( t )\rfloor$|> 0) in a one-year timespan (t ∈ [0,365] days). We note that the number of initial tumor cells is a constant parameter in our model and is independent of any other tumors that may be present in the individual (as we only consider the potential establishment of a new metastasis). Figure 3 shows examples of both successful metastatic tumors (Fig. 3A and C) and unsuccessful metastatic tumors (Fig. 3B and D) for different values of the MDSC delay τ (see Supplementary Text Section S7 for further description). For more examples of successful and unsuccessful tumors see Supplementary Figs. S7 and S8, respectively. While a metastatic tumor can become unsuccessful at any time point (and all tumors will be unsuccessful almost surely as t → ∞), we found that the tumor population is most likely to drop below one near the beginning of the simulation (i.e., soon after metastatic tumor seeding) when the tumor population is small (Supplementary Figs. S9 and S10).

Analysis of the probability of metastasis under different assumptions of MDSC–tumor–immune interactions for thousands of tumors studied in silico revealed striking dependencies of tumor outcomes on MDSC dynamics (Fig. 4). Through joint analysis of the effects of the number of circulating MDSCs (α₂) and the size of the MDSC delay (τ), we found that the probability of successful metastatic tumor establishment and the average size of metastatic tumors were positively correlated with the level of circulating MDSCs, and negatively correlated with the size of the MDSC delay. As more MDSCs became available at or near the site of the nascent metastasis, the NK-cell and CTL populations became more suppressed, resulting in a greater likelihood of tumor growth (Fig. 4A and B). The data indicate that the positive feedback loop (tumor cells are able to activate more MDSCs) reinforces the tumor's ability to grow, even in a “hot” tumor microenvironment.

We found that our model provides a biologically-driven way to determine exactly what can be inferred from levels of circulating MDSCs. Given the relative difficulty of defining MDSCs and the relative ease of sampling circulating cells this has clinical relevance (19). We saw that if the baseline level of circulating MDSCs (α₂) was high, MDSC activation delays had little effect on the metastasis establishment probability (Fig. 4A and B), but the MDSC delay still had a pronounced effect on the resulting sizes of the metastases that grew (Fig. 4C and D; Supplementary Fig. S11). Recall that we used a liberal definition for a successful metastasis: a population of > 1 tumor cells that survives for a year. Differences in the sizes of these nascent metastases from tens to thousands of cells bear direct clinical relevance. Further statistics on metastatic survival and size can be found in Supplementary Table S1. Relative to a MDSC delay of 0 days, a MDSC delay of 365 days led to a 2-fold decrease in the probability of successful metastasis, a 21-fold decrease in the mean tumor size (of successful tumors), and a 4.6-fold increase in the mean time to extinction of unsuccessful metastases (Fig. 4E and F).

Supplementary Figure S12 shows the effect of the rate of MDSC inhibition of NK cells (β₃) and Supplementary Fig. S13 shows the effect of the rate of MDSC inhibition of CTL cells (β₄). In these analyses, we found that more effective MDSCs (i.e., more inhibitory against antitumor populations) meant NK-cell and CTL populations were more inhibited, resulting in a larger tumor burden. However, if the level of inhibition of NK cells (β₃) was high enough, delays in recruitment of more MDSCs (τ) had little effect on the probability of successful metastatic tumors (as the tumor population will grow to very large levels very quickly, independently of a large increase in the number of MDSCs) but still affected the average size, as fewer NK cells resulted in more tumor cells (Supplementary Fig. S12). Because there are initially zero CTL cells and the CTL population does not reach extremely high levels relative to other populations (see for instance Fig. 1B–E, blue lines) changing β₄ did not have a large effect on the probability of successful metastasis (Supplementary Fig. S13A and S13B). However, increasing β₄ resulted in a small increase in the average size of successful tumors (see Supplementary Fig. S13C and S13D). We contrast this with analysis of CTL-abundant environments below in the Section titled “Analysis of tumor—MDSC dynamics in a CTL-abundant environment”.

MDSCs can be subdivided into one of two states: monocytic MDSCs (M-MDSCs; typically assumed to be more immunosuppressive) and granulocytic/poly-mononuclear (G- or PMN-MDSCs; refs. 1, 3, 6). The relative proportion of G- to M-MDSCs can alter the immunosuppressive properties of the tumor microenvironment (15, 73). For example, if the relative proportion of G- to M-MDSCs skews toward M-MDSCs, we would expect larger effects of MDSC delays (as seen in Fig. 4), whereas the opposite would be expected if G-MDSCs dominate. Extensions of the current model include separating M-MDSCs and G-MDSCs, with for instance β₃^M-MDSCs> β₃^G-MDSCs and β₄^M-MDSCs> β₄^G-MDSCs (see the Discussion for further details).

To summarize the results of this section, we have identified two crucial effects of MDSC delays on the stochastic tumor dynamics. First, that MDSC delays always resulted in significantly smaller tumor sizes. This effect became more pronounced as MDSCs became more immunesuppressive (in the model: when β₃, β₄ were large). Under these conditions, the increase in MDSCs allowed the tumor to outcompete the antitumor populations and reach large sizes. However if the MDSCs were so powerful as to completely inhibit the NK-cell and CTL populations, then increasing β₃, β₄ had no further effect. The effect of MDSC delay on tumor size was less pronounced when the MDSCs were less immunesuppressive (i.e., when β₃, β₄ are small): in this case increases in the number of MDSCs did not have significant effects on the long term dynamics of the other populations.

Second, that MDSC delays could result in drastically decreased probabilities of a successful new metastasis. This effect was most pronounced when the initial level of circulating MDSCs (α₂) was not too high, and when the MDSCs were not too immunesuppressive of the NK-cell population (large β₃). This was due to the greater likelihood of extinction of stochastic tumors (⁠|$\lfloor {x}_T( t )\rfloor$|< 1) early in the simulation. If the level of circulating MDSCs (α₂) was high, offering the nascent tumor protection against CTL and NK-cell responses, then the effects of delays in recruitment of more MDSCs were lessened. Similarly, if the MDSCs were strongly immunesuppressive (particularly against NK cells), then the tumor was likely to grow to a large size quickly, negating the impact of delays in MDSC recruitment on the probability of successful establishment of a new metastasis. These results establish how MDSC plasticity, as defined by their different suppressive functions and environments (i.e., circulation throughout the body or within a tumor), differentially contribute to tumor growth and progression of disease from a primary tumor location to a distant metastatic site.

To assess more rigorously the variability and uncertainty with which we know model parameters, we performed Bayesian parameter inference using clinical data on tumor progression as defined through RECIST (38). We fit the model, initialized in a lung microenvironment (Fig. 5A), to data from six individual tumors that spanned possible response outcomes (Fig. 5B). We selected a three-dimensional free parameter space (see Methods), consisting of the tumor growth rate, the NK-cell inhibition rate by MDSCs, and the rate of CTL stimulation by tumor–NK cell interactions. Successful fits were obtained for each of the tumors fit (Fig. 5C and D; Supplementary Figs. S14 and S15).

To analyze the parameters that give rise to different response dynamics, we plotted parameters sampled from the posteriors of each tumor fit (Fig. 5E-G). We saw a clear trend toward larger values of tumor growth rate (α₁) and NK-cell inhibition rate by MDSCs (β₃) for tumors that did not respond to treatment (tumors 5 and 6) compared with those that did respond to treatment (tumors 1 and 2; Fig. 5E). This can be understood in light of the previously characterized effects these parameters have on tumor growth (see Fig. 2D for example). Furthermore, strong correlations were observed for these parameters. The correlation between the two parameters was steeper for increasing tumors, suggestive of the discriminative ability of this parameter pair for quantifying tumor outcomes (i.e., whether tumors will grow or decay upon the initiation of treatment). In comparison, no correlations nor distinct effects on tumor outcomes were observed for the other two parameter pairs (Fig. 5F and G).

We tested the discriminative power of different combinations of posterior parameters by training decision trees to classify tumor responses as either decreasing (i.e., tumors 1 and 2) or increasing (i.e., tumors 5 and 6) over time. Supplementary Table S2 gives the cross validation scores for decision trees (maximum depth three) trained on different sets of posterior parameters as features. In line with the joint marginal posterior distributions observed (Fig. 5E) we saw that the best discriminative power was obtained using both the tumor growth rate (α₁) and the NK-cell inhibition rate by MDSCs (β₃) as features (81.5% mean accuracy). When constrained to using one feature, the NK-cell inhibition rate by MDSCs (69.1% mean accuracy) was a better predictor than the tumor growth rate (62.8% mean accuracy), even though the tumor growth rate was intricately tied to the classification outcome (43). Interest in interactions between MDSCs and NK cells has been growing in recent years (20, 21, 24); our results strongly suggest that more investigation is warranted.

Given that the clinical data did not capture immune dynamics, and the relative simplicity of the tumor dynamics in response to treatment, we expected to obtain fits to different individual tumor outcomes with the model. However, the relative importance of parameters that the fits revealed were unexpected. The strength of immunesuppressiveness—as controlled by NK-cell inhibition by MDSCs—was identified as the most important parameter in determining outcome. This has direct clinical implications: while it may not yet be possible to directly modulate this parameter in a clinical setting, it highlights the importance of interventions targeting properties of MDSCs in and around the tumor site. Moreover, successful fitting of various tumor responses to tumor—MDSC dynamics and the stratification of rate parameters that resulted demonstrates our ability to build and fit patient-specific tumor growth models (74), with which to predict metastatic outcomes.

To complement the analysis thus far of tumor–MDSC dynamics in NK cell–abundant environments, we studied how outcomes change if we considered an environment with relatively few NK cells but which was abundant in CTLs (see parameters in Methods and Fig. 6A). Simulation of the DDE model for different MDSC delays in the CTL-abundant environment (Supplementary Fig. S16) showed a dramatic increase in the steady state value of the CTL population, and declines in the NK population, relative to the NK cell–abundant environment (Fig. 1B–E). Because CTL cells must be activated before they act against the tumor Eq. (A), and since simulations start from the tumor-free steady state (initially zero CTL cells Eq. (F), then in the CTL-abundant environment small changes in the MDSC delay had negligible effects on the probability of establishing a new tumor (Fig. 6B and C; Supplementary Fig. S17). This highlights the importance of the environment: whereas small delays in MDSC recruitment/activation will impact tumor growth in the case that there are already circulating immune cells at the tumor site, this is not the case when the dominant immune component must be activated to respond to tumor growth.

We went on to fit the model to clinical response data, similar to above, but not assuming a CTL-abundant environment and adjusting the parameters accordingly. We observed an increase in the influence of α₆ (CTL stimulation by tumor–NK cell interaction) in classifying tumor outcomes (Fig. 6D–F; Supplementary Table S3; see also Supplementary Figs. S18 and S19 for model fits), thus CTL-related parameters became more important in this context (although note that the results of parameter inference in this setting are not directly comparable with those above due to differences in the prior distributions). The relative increase in importance makes sense intuitively: the increase in the CTL population size in a CTL-abundant environment can be expected to lead to greater prominence of this population in regulating tumor outcomes, and thus the predictive power of its parameters for classification. We observed a corresponding drop in the importance of NK cell–related parameters: the tumor growth rate was now the most important single parameter in determining tumor outcome (Supplementary Table S3); whereas previously the NK-cell inhibition rate by MDSCs was ranked highest.

Cancer dynamics are complex, and understanding cancer–immune dynamics is a complex systems biology problem (10, 28, 29, 39, 43). Modeling how tumors interact with the immune system is critical for understanding treatment responses and predicting the best possible therapeutic strategies in response to metastasis. MDSCs have been identified in various tumor microenvironments (8, 9, 15), where they can exert strong immunosuppressive effects leading to worse outcomes (12, 17, 37), yet a rigorous theoretical characterization of MDSC dynamics in the tumor microenvironment has remained lacking. Here, through the introduction of a SDDE model with which to study tumor–MDSC dynamics, we have provided means to characterize the plasticity of MDSCs and their effects on tumor progression and outcome.

With this model we began by studying outcomes under simple, idealized circumstances, such as: how large do tumors grow in the presence of MDSCs? What is their likelihood of persistence in the stochastic case? We discovered that delays in MDSC recruitment/activation had striking effects on metastatic growth and establishment. Under certain conditions (lower levels of circulating MDSCs), strategies that block MDSC recruitment to the site of the tumor are likely to greatly improve metastatic outcomes and hinder growth. We also demonstrated through model analyses how strategies that decrease the immunosuppressive properties of MDSCs are likely to have dramatic antitumor effects. Via Bayesian parameter estimation using data from tumor growth in vivo, we identified correlations between the tumor and the MDSC response parameters, again demonstrating the potential of inhibition of MDSCs as a desirable drug target.

Our inference results showed that in NK cell–abundant environments such as the lung, MDSC inhibition of NK cells was a crucial parameter determining outcomes; more important even than the tumor growth rate. We also found important differences between tumor microenvironments: we focused primarily on MDSC dynamics in the lung, an NK cell–rich environment (24, 46); however, we contrasted this with analysis of a CTL-abundant environments, that is, one in which CTL cells are greater in number than NK cells. In such a CTL-abundant environment, we observed a rise in prominence of the role of CTL activation (75, 76). Future work, where prompted by data, may need to examine intermediate environments, where CTL and NK-cell population sizes are more closely matched. Here, we might expect to see both NK-cell and CTL parameters to be important in informing tumor growth outcomes.

These results suggest that the identification of effective anti-MDSC treatment strategies to control cancer growth and spread ought to be more highly prioritized (8, 13, 17, 24). In particular, drug treatments that block MDSC recruitment to tumor sites and/or target MDSCs in the lymphoid organs might be expected to be most highly effective in preventing metastasis, but their effects are likely to be lessened if the level of circulating MDSCs is low or if MDSCs are less effective at suppressing antitumor populations. Because the level of circulating MDSCs (as well as the level of MDSC-immunosuppression) is likely to be highly variable within patients (20, 77), effective treatment strategies ought to be informed by patient-specific biomarkers (74, 78). In addition, evaluation of the phenotype of circulating MDSCs may not fully reflect the immunosuppressed state within tumors enough to predict potential response to immunotherapy, which may be determined in part by further mathematical and data-driven modeling. Toward this end, we have shown via tumor-specific parameter inference that we can train machine learning models using posterior parameters to classify metastatic outcomes. Future work, informed by more data (such as richer dynamic information or single-cell gene expression data) will provide additional means to classify treatment outcomes. In this context, it will be important to consider the prediction of responses in different tumor microenvironments and under different treatment regimes.

MDSCs cannot be assumed to be a homogeneous population. Although we have assumed as such here—for lack of data with which to quantify subpopulation-specific MDSC rate parameters—future models ought to consider MDSC heterogeneity. MDSCs are typically classified into one of two possible cell types, M-MDSCs and G-/PMN-MDSCs, which exhibit different levels of immunosuppression (1, 3, 15). M-MDSCs in metastatic breast cancer patients resemble monocytes isolated from patients with sepsis, indicating fascinating similarities between the immunosuppression capability of the MDSCs present in patients with metastatic (but perhaps not primary) breast cancer and those involved in the immunosuppressive sepsis response (79). Further measurement of MDSC subtype–specific immunosuppression in vivo will likely yield substantial new insight into their activity. Moreover, these additional data will permit the fitting of more detailed mathematical models that are able to describe patient-specific (or even tumor site-specific) dynamics and quantify the possible benefits of treatments targeting MDSCs. Current knowledge suggests that shifting MDSC phenotypes toward G-MDSCs is beneficial as this state is less immunosuppressive (1, 15); however, further characterization of these states is needed.

The models we have developed of MDSCs in the tumor microenvironment do not consider space, although of course spatial architectures play an important role in tumor progression (28, 32), in primary growth as well as for circulating tumor cells that seed metastases (34, 80). The role of spatial aspects of cancer niches in regulating MDSC–tumor dynamics will be an important topic in future work (81). Here, carefully fitting models to appropriate data ought to include both single cell–resolved characterization of the tumor microenvironment (15) and explicit spatial characterizations (82, 83). Furthermore, nonimmunologic mechanisms of MDSC suppression (e.g., tissue remodeling, angiogenesis, priming the metastatic niche (1)) may be important to consider in future models. Given our preclinical and clinical results (15), one current avenue of investigation is defining which mechanisms of MDSC suppression are affected by the epigenetic modulating drug entinostat, which sensitizes the lung metastatic tumor microenvironment and as such considers both immune- and nonimmune-mediated mechanisms of suppression. In addition, the level of suppression imparted by MDSCs within a tumor is not only directly affected by the signaling cascades and subsequent mechanism of suppression engaged (i.e., metabolic via amino acid depletion vs. immunologic suppression by IL10), but also indirectly affected by the surrounding microenvironments in which immune cell infiltration greatly varies. These are questions we hope to address in future work.

There is an urgent need to understand the role of MDSC dynamics during tumor growth and metastasis. Here, we discovered an essential and remarkable role for MDSC recruitment/activation in dictating growth outcomes in the context of new metastases. This is but the first step. To make progress, further conceptual model development tightly linked to inference and the gathering of higher-resolution data on MDSC phenotypes in vivo will be crucial. Mathematical modeling will continue to play an integral part in discovery as it allows us to account for the numerous and dynamic factors controlling MDSC plasticity and its impact of tumor responses in a way that traditional biologic biomarkers alone cannot. Only by developing theory and gathering data hand-in-hand can we hope to gain an understanding of the dynamics of MDSCs in the tumor microenvironment, and in turn, develop new therapies for metastatic disease.

E.T. Roussos Torres reports nonfinancial support from Synaptical outside the submitted work. A.L. MacLean reports grants from NIH NIGMS and grants from NSF DMS during the conduct of the study. No disclosures were reported by the other authors.

J. Kreger: Conceptualization, investigation, methodology, writing–original draft, writing–review and editing. E.T. Roussos Torres: Conceptualization, investigation, methodology, writing–review and editing. A.L. MacLean: Conceptualization, software, investigation, methodology, writing–original draft, writing–review and editing.

E.T. Roussos Torres acknowledges support from the Tower Cancer Research Foundation (Career Development Award), the Concern Foundation (Conquer Cancer Award), and the USC-Norris Comprehensive Cancer Center. A.L. MacLean acknowledges support from the NIH (R35GM143019) and the National Science Foundation (DMS2045327).

We would like to thank E.J. Fertig for valuable discussions and guidance, and the Tumor Microenvironment Program at the USC-Norris Comprehensive Cancer Center for their support. We would like to thank all members of the Roussos Torres and MacLean labs for valuable input and discussions. Figures 1A, 5A, and 6A were created with BioRender (BioRender.com).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.