Myeloid-derived suppressor cells (MDSC) play a prominent role in the tumor microenvironment. A quantitative understanding of the tumor–MDSC interactions that influence disease progression is critical, and currently lacking. We developed a mathematical model of metastatic growth and progression in immune-rich tumor microenvironments. We modeled the tumor–immune dynamics with stochastic delay differential equations and studied the impact of delays in MDSC activation/recruitment on tumor growth outcomes. In the lung environment, when the circulating level of MDSCs was low, the MDSC delay had a pronounced impact on the probability of new metastatic establishment: blocking MDSC recruitment could reduce the probability of metastasis by as much as 50%. To predict patient-specific MDSC responses, we fit to the model individual tumors treated with immune checkpoint inhibitors via Bayesian parameter inference. We reveal that control of the inhibition rate of natural killer (NK) cells by MDSCs had a larger influence on tumor outcomes than controlling the tumor growth rate directly. Posterior classification of tumor outcomes demonstrates that incorporating knowledge of the MDSC responses improved predictive accuracy from 63% to 82%. Investigation of the MDSC dynamics in an environment low in NK cells and abundant in cytotoxic T cells revealed, in contrast, that small MDSC delays no longer impacted metastatic growth dynamics. Our results illustrate the importance of MDSC dynamics in the tumor microenvironment overall and predict interventions promoting shifts toward less immune-suppressed states. We propose that there is a pressing need to consider MDSCs more often in analyses of tumor microenvironments.

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.

A stochastic delay differential equation model of tumor–immune dynamics in the presence of MDSCs

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.

Table 1.

Description of model parameters and values.

NotationDescriptionValueUnitsReferenceRange
xT(t),t ≤ 0 initial condition for tumor cells 1 or 2 — — — 
xMDSC(0) initial condition for MDSCs α22 — — — 
xNK(0) initial condition for NK cells |$\displaystyle \frac{{{\zeta }_2{\alpha }_4}}{{{\alpha }_2{\beta }_3 + {\zeta }_2{\zeta }_3}}$| — — — 
xCTL(0) initial condition for CTL cells — — — 
τ delay parameter for MDSCs varies days — — 
α1 tumor growth rate 10−1, varies days−1 (40, 84, 85[10−2,5 × 10−1
η tumor maximum size 107, varies — estimated [106,108
β1 tumor cells inhibition rate by NK cells 3.5 × 10−6 days−1 (40, 84, 85[10−7,10−6
β2 tumor cells inhibition rate by CTL cells 1.1 × 10−7 days−1 (84) [10−7,10−6
ζ1 tumor cell death rate 0, varies days−1 (30) [0,0.1] 
α2 MDSCs circulating rate 102, varies days−1 estimated (77[0,103
α3 MDSCs expansion coefficient 108 days−1 (23, 30, 86[107,109
ζ2 MDSCs death rate 0.2 days−1 (87, 88[0,1] 
α4 NK cells circulating rate 1.4 × 104, varies days−1 (84[103,105
α5 NK cells expansion coefficient 2.5 × 10−2 days−1 (40, 84, 85[10−2,10−1
β3 NK cells inhibition rate by MDSCs 4 × 10−5, varies days−1 (30[10−5,10−4
ζ3 NK cells death rate 4.12 × 10−2 days−1 (84[10−2,10−1
α6 CTL stimulation by tumor-NK cell interaction 1.1 × 10−7, varies days−1 (89, 90[10−7,10−6
α7 CTL expansion coefficient 10−1 days−1 (91[5 × 10−2,10−1
β4 CTL inhibition rate by MDSCs 10−4, varies days−1 (30[5 × 10−5,5 × 10−4
ζ4 CTL death rate 2 × 10−2 days−1 (40, 89[10−2,10−1
γ1 steepness of MDSC production 1010 — (30, 86[109,1011
γ2 steepness of NK production 2.02 × 107 — (40, 84[106,108
γ3 steepness of CTL production 2.02 × 107 — (40, 84, 85[106,108
NotationDescriptionValueUnitsReferenceRange
xT(t),t ≤ 0 initial condition for tumor cells 1 or 2 — — — 
xMDSC(0) initial condition for MDSCs α22 — — — 
xNK(0) initial condition for NK cells |$\displaystyle \frac{{{\zeta }_2{\alpha }_4}}{{{\alpha }_2{\beta }_3 + {\zeta }_2{\zeta }_3}}$| — — — 
xCTL(0) initial condition for CTL cells — — — 
τ delay parameter for MDSCs varies days — — 
α1 tumor growth rate 10−1, varies days−1 (40, 84, 85[10−2,5 × 10−1
η tumor maximum size 107, varies — estimated [106,108
β1 tumor cells inhibition rate by NK cells 3.5 × 10−6 days−1 (40, 84, 85[10−7,10−6
β2 tumor cells inhibition rate by CTL cells 1.1 × 10−7 days−1 (84) [10−7,10−6
ζ1 tumor cell death rate 0, varies days−1 (30) [0,0.1] 
α2 MDSCs circulating rate 102, varies days−1 estimated (77[0,103
α3 MDSCs expansion coefficient 108 days−1 (23, 30, 86[107,109
ζ2 MDSCs death rate 0.2 days−1 (87, 88[0,1] 
α4 NK cells circulating rate 1.4 × 104, varies days−1 (84[103,105
α5 NK cells expansion coefficient 2.5 × 10−2 days−1 (40, 84, 85[10−2,10−1
β3 NK cells inhibition rate by MDSCs 4 × 10−5, varies days−1 (30[10−5,10−4
ζ3 NK cells death rate 4.12 × 10−2 days−1 (84[10−2,10−1
α6 CTL stimulation by tumor-NK cell interaction 1.1 × 10−7, varies days−1 (89, 90[10−7,10−6
α7 CTL expansion coefficient 10−1 days−1 (91[5 × 10−2,10−1
β4 CTL inhibition rate by MDSCs 10−4, varies days−1 (30[5 × 10−5,5 × 10−4
ζ4 CTL death rate 2 × 10−2 days−1 (40, 89[10−2,10−1
γ1 steepness of MDSC production 1010 — (30, 86[109,1011
γ2 steepness of NK production 2.02 × 107 — (40, 84[106,108
γ3 steepness of CTL production 2.02 × 107 — (40, 84, 85[106,108

Note: Estimated from the literature (see in particular refs. 40, 77, 84, 86). Cell populations are measured in terms of cell numbers and are nondimensionalized. The first column is the parameter notation, the second column is the parameter description, the third column is the parameter estimated value, the fourth column is the parameter units (if applicable), the fifth column is the citation of the reference for the parameter estimate, and the sixth column is the parameter range used for the global stability analysis (see Materials and Methods).

Figure 1.

Larger MDSC delays result in significantly altered tumor growth dynamics. A, Model diagram. The MDSC, NK cell, and CTL populations are all signaled to proliferate in the presence of a metastatic tumor. The MDSC population inhibits the NK-cell and CTL populations, and the NK-cell and CTL populations inhibit the tumor population. B–D, Simulations of the deterministic (DDE) system [Eq. (C), g = 0] over one year, with one initial tumor cell and different MDSC delay parameter τ. See Materials and Methods for simulation details. B,τ = 0. C,τ = 10. D,τ = 50. E,τ = 365.

Figure 1.

Larger MDSC delays result in significantly altered tumor growth dynamics. A, Model diagram. The MDSC, NK cell, and CTL populations are all signaled to proliferate in the presence of a metastatic tumor. The MDSC population inhibits the NK-cell and CTL populations, and the NK-cell and CTL populations inhibit the tumor population. B–D, Simulations of the deterministic (DDE) system [Eq. (C), g = 0] over one year, with one initial tumor cell and different MDSC delay parameter τ. See Materials and Methods for simulation details. B,τ = 0. C,τ = 10. D,τ = 50. E,τ = 365.

Close modal

We denote xT, xMDSC, xNK, and xCTL 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 δxi denotes the rate of change of xi, i ∈ [T, MDSC, NK, CTL].

formula
On the basis of these biological processes, we developed a stochastic delay differential equation (SDDE) model to characterize tumor–immune interactions. This model is described by the equation:
formula
at time t, with delay 0 < τ < t, where f(·) describes the deterministic dynamics controlled by the model interactions, g(·)dW(t) describes the stochastic dynamics, dW(t) denotes an increment of a Weiner process, W(t), and xj(t) = [xT(t),xMDSC(t),xNK(t),xCTL(t)]. The model thus consists of coupled SDDEs, where we assumed an Itô interpretation (49). For the stochastic dynamics, we have:
formula

where ξi(t) is the size of the ith 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:

formula

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 β1 and β2, respectively. MDSCs are activated due to their basal circulation, α2, and die at rate ζ2. In the presence of tumor cells, immune-suppressive signals (α3) 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, α3 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 xMDSC; 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 α4, and CTLs are activated by the NK cell–tumor cell interaction at rate α6. In line with (30), both NK cells and CTLs can be activated by the tumor (at rates α5 and α7, respectively). We assumed that NK cells and CTLs can be inhibited by MDSCs (at rates β3 and β4, respectively), and are lost due to cell death (at rates ζ3 and ζ4, 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

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.

Bayesian parameter inference with RECIST data

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 (mm3) as |$\frac{1}{2}{x}^3\$|(38, 64). We estimated the number of tumor cells from this volume by multiplying by a factor of 107 (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(xT +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: β3 (NK cells inhibition rate by MDSCs), α6 (CTL stimulation by tumor–NK cell interaction), and α1 (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 (xT(0) = 8,395.4, xMDSC(0) = 804.1, xNK(0) = 197,565.7 and xCTL(0) = 1,654.4). Therefore, we rescaled η = 105 (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:

formula

where σ is the observational noise. For each tumor, we ran four independent Markov chain Monte Carlo (MCMC) simulations with 2 × 103 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 xT(0) = 16,654.8, xMDSC(0) = 707.8, xNK(0) = 16,162.1 and xCTL(0) = 778,128.8 (simulation values at day ten) and that the prior distributions for the parameters were:

formula

Decision tree classification of tumor responses

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).

Modeling MDSC dynamics in a CTL-abundant environment

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 α4 = 103, and increased the CTL stimulation by tumor/NK cell–interaction rate to α6 = 2 × 10−3. This led to an increase in the number of CTL cells and a concordant decrease in the number of NK cells.

Data availability statement

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.

Dynamics of metastatic growth in the presence of MDSCs

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. 1BE). 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 (xT = 0), the tumor-free fixed point of the model is:

formula
where |$\hat{x}$|T, |$\hat{x}$|MDSC, |$\hat{x}$|NK, and |$\hat{x}$|CTL represent the steady state values of xT, xMDSC, xNK, and xCTL, respectively. We observed baseline populations of MDSCs and NK cells at the metastatic site, but no CTL cells, as they need to be recruited and activated against the tumor. Because tumor cells cannot be spontaneously generated in this model, the tumor-free fixed point Eq. (F) was stable (Supplementary Text Section S3). In the case of a nonzero tumor population ( |$\hat{x}$|T> 0), in general the steady state must be determined numerically, although we can derive analytic approximations in special cases. For example, for |$\hat{x}$|T> 0, the steady states of the nontumor populations are:
formula
where
formula

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 xT> 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 R0 in epidemiological models; see Supplementary Text Section S5 and Supplementary Fig. S5 for details). G is given by:

formula

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×106. We changed the tumor cell death rate (ζ1) 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).

Figure 2.

Dependencies of tumor growth characteristics on model parameters. Simulations of the ODE system [Eq. (C), g = τ = 0] with one initial tumor cell. A–C, Different tumor growth thresholds G Eq. (H). A, G ≈ 0.8 (parameters as in Table 1). B, G ≈ 0 (ζ1 = 0.81). C, G ≈ −0.2 (ζ1 = 1). D, Morris global sensitivity analysis (GSA) for the steady state of the tumor population for all model parameters. Green denotes parameters that are positively correlated with the tumor size at steady state; red denotes negatively correlated. Hexagons represent MDSC-specific parameters; circles represent non–MDSC-specific parameters. E and F, Effects of the NK-cell inhibition rate by MDSCs (β3), for β3 = 10−5, the minimum of the GSA range (E); the tumor size at steady state is 2.5 × 102. And for β3 = 10−4, the maximum of the GSA range (F); The tumor size at steady state is 9.9 × 106.

Figure 2.

Dependencies of tumor growth characteristics on model parameters. Simulations of the ODE system [Eq. (C), g = τ = 0] with one initial tumor cell. A–C, Different tumor growth thresholds G Eq. (H). A, G ≈ 0.8 (parameters as in Table 1). B, G ≈ 0 (ζ1 = 0.81). C, G ≈ −0.2 (ζ1 = 1). D, Morris global sensitivity analysis (GSA) for the steady state of the tumor population for all model parameters. Green denotes parameters that are positively correlated with the tumor size at steady state; red denotes negatively correlated. Hexagons represent MDSC-specific parameters; circles represent non–MDSC-specific parameters. E and F, Effects of the NK-cell inhibition rate by MDSCs (β3), for β3 = 10−5, the minimum of the GSA range (E); the tumor size at steady state is 2.5 × 102. And for β3 = 10−4, the maximum of the GSA range (F); The tumor size at steady state is 9.9 × 106.

Close modal

Parameter sensitivity analysis reveals that inhibition rates between populations are important in determining tumor growth outcomes

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 α2, α3, ζ2, β3, β4, and γ1. 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 ζ2 (death rate of MDSCs) and γ1 (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, β3 (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. 1BE and Eq. (F)]: the MDSC population must effectively suppress NK cells for the tumor to be able to grow rather than die out. Similarly, β4 (inhibition of CTL cells by MDSCs) was important, but less so than β3 as the CTL population is initially small thus less important for initial growth of the tumor [see CTL cells in Fig. 1BE 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 β3 (inhibition of NK cells by MDSCs) on the tumor steady state at both ends of the range studied. We saw that small values of β3 (Fig. 2E) resulted in small tumors (tumor steady state approx. 2.5×102), whereas large values of β3 (Fig. 2F) resulted in much larger tumors (tumor steady state approx. 9.9 × 106).

Comparison of MDSC-specific parameters (hexagons) with non–MDSC-specific parameters (circles, Fig. 2D) showed that α2, α3, α1, η, β3, ζ3, β4, ζ4, γ2, and γ3 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 α6, β1, β2, β3, and β4. α6 is the rate of CTL stimulation by tumor–NK cell interaction, β1 and β2 are inhibition rates of tumor cells by NK and CTL cells, and β3 and β4 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).

Stochastic dynamics of metastatic growth and establishment

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).

Figure 3.

Stochastic effects influence the growth and probability of establishment of metastatic tumors. Examples of simulations of the SDDE system Eq. (C) over one year, with two initial tumor cells and different values of the MDSC delay parameter, τ (all other parameters not explicitly defined otherwise are as in Table 1). A “successful” metastatic tumor is one that does not drop below a size of one tumor cell over the simulation period. A,τ = 0; successful. B,τ = 10; unsuccessful. C,τ = 50; successful. D,τ = 365; unsuccessful.

Figure 3.

Stochastic effects influence the growth and probability of establishment of metastatic tumors. Examples of simulations of the SDDE system Eq. (C) over one year, with two initial tumor cells and different values of the MDSC delay parameter, τ (all other parameters not explicitly defined otherwise are as in Table 1). A “successful” metastatic tumor is one that does not drop below a size of one tumor cell over the simulation period. A,τ = 0; successful. B,τ = 10; unsuccessful. C,τ = 50; successful. D,τ = 365; unsuccessful.

Close modal

Delays in MDSC recruitment decrease the probability of metastasis and the size of metastatic tumors

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 (α2) 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.

Figure 4.

Effects of MDSC properties on the probability of metastatic establishment. Stochastic simulations run for a period of one year. Each point is the mean over at least 105 simulations. All parameters not explicitly defined otherwise are as in Table 1. Ribbons (shaded area) represent the SE. A, Probability of new tumor establishment over a period of one year, for different values of the level of circulating MDSCs (α2) and the MDSC delay (τ). B, As for A with τ plotted on log scale. C, Of the new metastases that are successfully established, the distribution of their means sizes is given. D, As for C with τ plotted on log scale. E, Of the new metastases that go extinct, the distribution of the mean times to extinction is given. F, As for E with τ plotted on log scale.

Figure 4.

Effects of MDSC properties on the probability of metastatic establishment. Stochastic simulations run for a period of one year. Each point is the mean over at least 105 simulations. All parameters not explicitly defined otherwise are as in Table 1. Ribbons (shaded area) represent the SE. A, Probability of new tumor establishment over a period of one year, for different values of the level of circulating MDSCs (α2) and the MDSC delay (τ). B, As for A with τ plotted on log scale. C, Of the new metastases that are successfully established, the distribution of their means sizes is given. D, As for C with τ plotted on log scale. E, Of the new metastases that go extinct, the distribution of the mean times to extinction is given. F, As for E with τ plotted on log scale.

Close modal

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 (α2) 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 (β3) and Supplementary Fig. S13 shows the effect of the rate of MDSC inhibition of CTL cells (β4). 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 (β3) 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. 1BE, blue lines) changing β4 did not have a large effect on the probability of successful metastasis (Supplementary Fig. S13A and S13B). However, increasing β4 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 β3M-MDSCs> β3G-MDSCs and β4M-MDSCs> β4G-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 β3, β4 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 β3, β4 had no further effect. The effect of MDSC delay on tumor size was less pronounced when the MDSCs were less immunesuppressive (i.e., when β3, β4 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 (α2) was not too high, and when the MDSCs were not too immunesuppressive of the NK-cell population (large β3). 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 (α2) 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.

Bayesian parameter inference reveals the clinical importance of MDSC–NK cell interactions in the lung microenvironment

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).

Figure 5.

Interactions between MDSCs and NK cells control clinical tumor growth outcomes. All parameters not explicitly defined otherwise are as in Table 1. A, Schematic diagram of NK cell–abundant environment. B, Relative change in tumor size from the baseline assessment for six tumors from patients with non–small cell lung cancer undergoing treatment with anti–PD-L1. Tumors are ordered (1–6) by their response, compared to baseline assessment. C, Tumor 2 model trajectories based on the relative change in the tumor population with the black dots representing the data, the purple line representing the fit from using the median of the posterior distribution for each parameter, and the shaded area denoting the 90% credible interval (where 90% of the posterior trajectories lie). D, Same as B for tumor 5. E–G, Samples from the posterior distribution of each of the six tumors, 8×103 samples plotted for pairs of model parameters. E, NK-cell inhibition rate by MDSCs (β3) versus tumor growth rate (α1). F, NK-cell inhibition rate by MDSCs (β3) versus CTL stimulation by tumor–NK cell interaction (α6) G: CTL stimulation by tumor–NK cell interaction (α6) versus tumor growth rate (α1).

Figure 5.

Interactions between MDSCs and NK cells control clinical tumor growth outcomes. All parameters not explicitly defined otherwise are as in Table 1. A, Schematic diagram of NK cell–abundant environment. B, Relative change in tumor size from the baseline assessment for six tumors from patients with non–small cell lung cancer undergoing treatment with anti–PD-L1. Tumors are ordered (1–6) by their response, compared to baseline assessment. C, Tumor 2 model trajectories based on the relative change in the tumor population with the black dots representing the data, the purple line representing the fit from using the median of the posterior distribution for each parameter, and the shaded area denoting the 90% credible interval (where 90% of the posterior trajectories lie). D, Same as B for tumor 5. E–G, Samples from the posterior distribution of each of the six tumors, 8×103 samples plotted for pairs of model parameters. E, NK-cell inhibition rate by MDSCs (β3) versus tumor growth rate (α1). F, NK-cell inhibition rate by MDSCs (β3) versus CTL stimulation by tumor–NK cell interaction (α6) G: CTL stimulation by tumor–NK cell interaction (α6) versus tumor growth rate (α1).

Close modal

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 (α1) and NK-cell inhibition rate by MDSCs (β3) 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 (α1) and the NK-cell inhibition rate by MDSCs (β3) 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.

Analysis of tumor—MDSC dynamics in a CTL-abundant environment

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. 1BE). 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.

Figure 6.

The impact of CTL parameters on tumor growth increases in a CTL-abundant environment. A, Schematic diagram of CTL-abundant environment. B, Stochastic simulations run for a period of one year. Each point is the mean over at least 105 simulations. Ribbons (shaded area) represent the SE. Probability of new tumor establishment over a period of one year, for different values of the level of circulating MDSCs (α2) and the MDSC delay (τ). C, As for A with τ ∈ [0,10]. Simulations use parameters as defined in the final section of the Methods, all other parameters not explicitly defined otherwise are as in Table 1. D–F, Samples from the posterior distribution of each of the six tumors, 8 × 103 samples plotted for pairs of model parameters. D, NK-cell inhibition rate by MDSCs (β3) versus tumor growth rate (α1). E, NK-cell inhibition rate by MDSCs (β3) versus CTL stimulation by tumor–NK cell interaction (α6). F, CTL stimulation by tumor–NK cell interaction (α6) versus tumor growth rate (α1).

Figure 6.

The impact of CTL parameters on tumor growth increases in a CTL-abundant environment. A, Schematic diagram of CTL-abundant environment. B, Stochastic simulations run for a period of one year. Each point is the mean over at least 105 simulations. Ribbons (shaded area) represent the SE. Probability of new tumor establishment over a period of one year, for different values of the level of circulating MDSCs (α2) and the MDSC delay (τ). C, As for A with τ ∈ [0,10]. Simulations use parameters as defined in the final section of the Methods, all other parameters not explicitly defined otherwise are as in Table 1. D–F, Samples from the posterior distribution of each of the six tumors, 8 × 103 samples plotted for pairs of model parameters. D, NK-cell inhibition rate by MDSCs (β3) versus tumor growth rate (α1). E, NK-cell inhibition rate by MDSCs (β3) versus CTL stimulation by tumor–NK cell interaction (α6). F, CTL stimulation by tumor–NK cell interaction (α6) versus tumor growth rate (α1).

Close modal

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 α6 (CTL stimulation by tumor–NK cell interaction) in classifying tumor outcomes (Fig. 6DF; 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.

Note: Supplementary data for this article are available at Cancer Immunology Research Online (http://cancerimmunolres.aacrjournals.org/).

1.
Veglia
F
,
Sanseviero
E
,
Gabrilovich
DI
.
Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity
.
Nat Rev Immunol
2021
;
21
:
485
98
.
2.
Hegde
S
,
Leader
AM
,
Merad
M
.
MDSC: Markers, development, states, and unaddressed complexity
.
Immunity
2021
;
54
:
875
84
.
3.
Bergenfelz
C
,
Leandersson
K
.
The generation and identity of human myeloid-derived suppressor cells
.
Front Oncol
2020
;
10
:
109
.
4.
Talmadge
JE
,
Gabrilovich
DI
.
History of myeloid-derived suppressor cells (MDSCs) in the macro- and micro-environment of tumor-bearing hosts
.
Nat Rev Cancer
2013
;
13
:
739
52
.
5.
Bronte
V
,
Brandau
S
,
Chen
S-H
,
Colombo
MP
,
Frey
AB
,
Greten
TF
, et al
.
Recommendations for myeloid-derived suppressor cell nomenclature and characterization standards
.
Nat Commun
2016
;
7
:
12150
.
6.
Veglia
F
,
Perego
M
,
Gabrilovich
D
.
Myeloid-derived suppressor cells coming of age
.
Nat Immunol
2018
;
19
:
108
19
.
7.
Alshetaiwi
H
,
Pervolarakis
N
,
McIntyre
LL
,
Ma
D
,
Nguyen
Q
,
Rath
JA
, et al
.
Defining the emergence of myeloid-derived suppressor cells in breast cancer using single-cell transcriptomics
.
Sci Immunol
2020
;
5
:
eaay6017
.
8.
Tesi
RJ
.
MDSC; the most important cell you have never heard of
.
Trends Pharmacol Sci
2019
;
40
:
4
7
.
9.
Ren
L
,
Li
J
,
Wang
C
,
Lou
Z
,
Gao
S
,
Zhao
L
, et al
.
Single cell RNA sequencing for breast cancer: present and future
.
Cell Death Discov
2021
;
7
:
104
.
10.
Hamilton
PT
,
Anholt
BR
,
Nelson
BH
.
Tumour immunotherapy: lessons from predator–prey theory
.
Nat Rev Immunol
2022
;
22
:
765
75
.
11.
Wynn
TA
.
Myeloid-cell differentiation redefined in cancer
.
Nat Immunol
2013
;
14
:
197
9
.
12.
Cha
YJ
,
Koo
JS
.
Role of tumor-associated myeloid cells in breast cancer
.
Cells
2020
;
9
:
1785
.
13.
Kumar
V
,
Patel
S
,
Tcyganov
E
,
Gabrilovich
DI
.
The nature of myeloid-derived suppressor cells in the tumor microenvironment
.
Trends Immunol
2016
;
37
:
208
20
.
14.
Tcyganov
E
,
Mastio
J
,
Chen
E
,
Gabrilovich
DI
.
Plasticity of myeloid-derived suppressor cells in cancer
.
Curr Opin Immunol
2018
;
51
:
76
82
.
15.
Sidiropoulos
DN
,
Rafie
CI
,
Jang
JK
,
Castanon
S
,
Baugh
AG
,
Gonzalez
E
, et al
.
Entinostat decreases immune suppression to promote antitumor responses in a HER2+ breast tumor microenvironment
.
Cancer Immunol Res
2022
;
10
:
656
69
.
16.
Christmas
BJ
,
Rafie
CI
,
Hopkins
AC
,
Scott
BA
,
Ma
HS
,
Cruz
KA
, et al
.
Entinostat converts immune-resistant breast and pancreatic cancers into checkpoint-responsive tumors by reprogramming tumor-infiltrating MDSCs
.
Cancer Immunol Res
2018
;
6
:
1561
77
.
17.
De Cicco
P
,
Ercolano
G
,
Ianaro
A
.
The new era of cancer immunotherapy: targeting myeloid-derived suppressor cells to overcome immune evasion
.
Front Immunol
2020
;
11
:
1680
.
18.
Gonda
K
,
Shibata
M
,
Ohtake
T
,
Matsumoto
Y
,
Tachibana
K
,
Abe
N
, et al
.
Myeloid-derived suppressor cells are increased and correlated with type 2 immune responses, malnutrition, inflammation, and poor prognosis in patients with breast cancer
.
Oncol Lett
2017
;
14
:
1766
74
.
19.
Angell
TE
,
Lechner
MG
,
Smith
AM
,
Martin
SE
,
Groshen
SG
,
Maceri
DR
, et al
.
Circulating myeloid-derived suppressor cells predict differentiated thyroid cancer diagnosis and extent
.
Thyroid
2016
;
26
:
381
9
.
20.
Tumino
N
,
Di Pace
AL
,
Besi
F
,
Quatrini
L
,
Vacca
P
,
Moretta
L
.
Interaction between MDSC and NK cells in solid and hematological malignancies: impact on HSCT
.
Front Immunol
2021
;
12
:
316
.
21.
Bruno
A
,
Mortara
L
,
Baci
D
,
Noonan
DM
,
Albini
A
.
Myeloid derived suppressor cells interactions with natural killer cells and pro-angiogenic activities: roles in tumor progression
.
Front Immunol
2019
;
10
:
771
.
22.
Hoechst
B
,
Voigtlaender
T
,
Ormandy
L
,
Gamrekelashvili
J
,
Zhao
F
,
Wedemeyer
H
, et al
.
Myeloid derived suppressor cells inhibit natural killer cells in patients with hepatocellular carcinoma via the NKp30 receptor
.
Hepatology
2009
;
50
:
799
807
.
23.
Cao
Y
,
Feng
Y
,
Zhang
Y
,
Zhu
X
,
Jin
F
.
L-Arginine supplementation inhibits the growth of breast cancer by enhancing innate and adaptive immune responses mediated by suppression of MDSCs in vivo
.
BMC Cancer
2016
;
16
:
343
.
24.
Zalfa
C
,
Paust
S
.
Natural killer cell interactions with myeloid derived suppressor cells in the tumor microenvironment and implications for cancer immunotherapy
.
Front Immunol
2021
;
12
:
633205
.
25.
Markowitz
J
,
Wesolowski
R
,
Papenfuss
T
,
Brooks
TR
,
Carson
WE
.
Myeloid-derived suppressor cells in breast cancer
.
Breast Cancer Res Treat
2013
;
140
:
13
21
.
26.
Jarrett
AM
,
Shah
A
,
Bloom
MJ
,
McKenna
MT
,
Hormuth
DA
,
Yankeelov
TE
, et al
.
Experimentally-driven mathematical modeling to improve combination targeted and cytotoxic therapy for HER2+ breast cancer
.
Sci Rep
2019
;
9
:
12830
.
27.
Lai
X
,
Stiff
A
,
Duggan
M
,
Wesolowski
R
,
Carson
WE
,
Friedman
A
.
Modeling combination therapy for breast cancer with BET and immune checkpoint inhibitors
.
Proc Nat Acad Sci USA
2018
;
115
:
5534
9
.
28.
Mahlbacher
GE
,
Reihmer
KC
,
Frieboes
HB
.
Mathematical modeling of tumorimmune cell interactions
.
J Theor Biol
2019
;
469
:
47
60
.
29.
Louzoun
Y
.
The evolution of mathematical immunology
.
Immunol Rev
2007
;
216
:
9
20
.
30.
Shariatpanahi
SP
,
Shariatpanahi
SP
,
Madjidzadeh
K
,
Hassan
M
,
Abedi-Valugerdi
M
.
Mathematical modeling of tumor-induced immunosuppression by myeloid-derived suppressor cells: implications for therapeutic targeting strategies
.
J Theor Biol
2018
;
442
:
1
10
.
31.
Allahverdy
A
,
Moghaddam
AK
,
Rahbar
S
,
Shafiekhani
S
,
Mirzaie
HR
,
Amanpour
S
, et al
.
An agent-based model for investigating the effect of myeloid-derived suppressor cells and its depletion on tumor immune surveillance
.
J Med Signals Sen
2019
;
9
:
15
.
32.
Liao
K-L
,
Bai
X-F
,
Friedman
A
.
Mathematical modeling of interleukin-35 promoting tumor growth and angiogenesis
.
PLoS One
2014
;
9
:
e110126
.
33.
Siewe
N
,
Friedman
A
.
TGF-β inhibition can overcome cancer primary resistance to PD-1 blockade: A mathematical model
.
PLoS One
2021
;
16
:
e0252620
.
34.
Mukherjee
M
,
Levine
H
.
Cluster size distribution of cells disseminating from a primary tumor
.
PLoS Comput Biol
2021
;
17
:
e1009011
.
35.
Lai
S-CA
,
Gundlapalli
H
,
Ekiz
HA
,
Jiang
A
,
Fernandez
E
,
Welm
AL
.
Blocking short-form ron eliminates breast cancer metastases through accumulation of stem-like CD4+ T cells that subvert immunosuppression
.
Cancer Discov
2021
;
11
:
3178
97
.
36.
Gui
P
,
Bivona
TG
.
Evolution of metastasis: new tools and insights
.
Trends Cancer
2022
;
8
:
98
109
.
37.
Trovato
R
,
Canè
S
,
Petrova
V
,
Sartoris
S
,
Ugel
S
,
De Sanctis
F
.
The engagement between MDSCs and metastases: partners in crime
.
Front Oncol
2020
;
10
:
165
.
38.
Ghaffari Laleh
N
,
Loeffler
CML
,
Grajek
J
,
Staňková
K
,
Pearson
AT
,
Muti
HS
, et al
.
Classical mathematical models for prediction of response to chemotherapy and immunotherapy
.
PLoS Comput Biol
2022
;
18
:
e1009822
.
39.
Wodarz
D
,
Komarova
N
.
Dynamics of cancer: mathematical foundations Of oncology.
World Scientific
;
2014
.
40.
de Pillis
LG
,
Radunskaya
AE
,
Wiseman
CL
.
A validated mathematical model of cell-mediated immune response to tumor growth
.
Cancer Res
2005
;
65
:
7950
8
.
41.
Rodriguez-Brenes
IA
,
Komarova
NL
,
Wodarz
D
.
Tumor growth dynamics: insights into evolutionary processes
.
Trends Ecol Evol
2013
;
28
:
597
604
.
42.
Zahir
N
,
Sun
R
,
Gallahan
D
,
Gatenby
RA
,
Curtis
C
.
Characterizing the ecological and evolutionary dynamics of cancer
.
Nat Genet
2020
;
52
:
759
67
.
43.
Sontag
ED
.
A dynamic model of immune responses to antigen presentation predicts different regions of tumor or pathogen elimination
.
Cell Syst
2017
;
4
:
231
41
.
44.
Sihto
H
,
Lundin
J
,
Lundin
M
,
Lehtimäki
T
,
Ristimäki
A
,
Holli
K
, et al
.
Breast cancer biological subtypes and protein expression predict for the preferential distant metastasis sites: a nationwide cohort study
.
Breast Cancer Res
2011
;
13
:
R87
.
45.
Xiao
W
,
Zheng
S
,
Liu
P
,
Zou
Y
,
Xie
X
,
Yu
P
, et al
.
Risk factors and survival outcomes in patients with breast cancer and lung metastasis: a population-based study
.
Cancer Med
2018
;
7
:
922
30
.
46.
Hervier
B
,
Russick
J
,
Cremer
I
,
Vieillard
V
.
NK cells in the human lungs
.
Front Immunol
2019
;
10
:
1263
.
47.
Chan
IS
,
Ewald
AJ
.
The changing role of natural killer cells in cancer metastasis
.
J Clin Invest
2022
;
132
:
e143762
.
48.
Chan
IS
,
Knútsdóttir
H
,
Ramakrishnan
G
,
Padmanaban
V
,
Warrier
M
,
Ramirez
JC
, et al
.
Cancer cells educate natural killer cells to a metastasis-promoting cell state
.
J Cell Biol
2020
;
219
:
e202001134
.
49.
Øksendal
B
.
Stochastic Differential Equations: An Introduction with Applications
. 6th ed. (
Springer
,
2013
).
50.
Fox
RF
.
Functional-calculus approach to stochastic differential equations
.
Phys Rev A
1986
;
33
:
467
76
.
51.
Coomer
MA
,
Ham
L
,
Stumpf
MPH
.
Noise distorts the epigenetic landscape and shapes cell-fate decisions
.
Cell Syst
2022
;
13
:
83
102
.
52.
Wu
C
,
Hua
Q
,
Zheng
L
.
Generation of myeloid cells in cancer: the spleen matters
.
Front Immunol
2020
;
11
:
1126
.
53.
Mackey
M
,
Glass
L
.
Oscillation and chaos in physiological control systems
.
Science
1977
;
197
:
287
9
.
54.
Bezanson
J
,
Edelman
A
,
Karpinski
S
,
Shah
VB
.
Julia: a fresh approach to numerical computing
.
SIAM Rev
2017
;
59
:
65
98
.
55.
Rackauckas
C
,
Nie
Q
.
Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia
.
J Open Res Softw
2017
;
5
.
56.
Rackauckas
C
,
Nie
Q
.
Stability-optimized high order methods and stiffness detection for pathwise stiff stochastic differential equations
.
In
:
2020 IEEE High Performance Extreme Computing Conference (HPEC)
,
1
8
.
Waltham, MA: IEEE;
2020
.
Available from:
https://ieeexplore.ieee.org/document/9286178/.
57.
Roesch
E
,
Greener
JG
,
MacLean
AL
,
Nassar
H
,
Rackauckas
C
,
Holy
TE
, et al
.
Julia for biologists
.
arXiv
:
2021
,
2109.09973 [q-bio] Available from:
http://arxiv.org/abs/2109.09973. ArXiv: 2109.09973.
58.
Saltelli
AA
.
Sensitivity analysis in practice: a guide to assessing scientific models
.
Chichester/Hoboken, NJ
:
Wiley
;
2004
.
Available from:
http://archive.org/details/sensitivityanaly00salt.
59.
Liu
D
,
Li
L
,
Rostami-Hodjegan
A
,
Bois
FY
,
Jamei
M
.
Considerations and caveats when applying global sensitivity analysis methods to physiologically based pharmacokinetic models
.
AAPS J
2020
;
22
:
93
.
60.
Therasse
P
,
Arbuck
SG
,
Eisenhauer
EA
,
Wanders
J
,
Kaplan
RS
,
Rubinstein
L
, et al
.
New guidelines to evaluate the response to treatment in solid tumors
.
J Natl Cancer Inst
2000
;
92
:
205
16
.
61.
Spigel
DR
,
Chaft
JE
,
Gettinger
S
,
Chao
BH
,
Dirix
L
,
Schmid
P
, et al
.
FIR: efficacy, safety, and biomarker analysis of a Phase II open-label study of atezolizumab in PD-L1–selected patients with NSCLC
.
J Thorac Oncol
2018
;
13
:
1733
42
.
62.
Linden
Wvd
,
Dose
V
,
Toussaint
Uv
.
Bayesian probability theory: applications in the physical sciences
.
Cambridge, United Kingdom:
Cambridge University Press
;
2014
.
Available from:
https://www.cambridge.org/core/books/bayesian-probability-theory/7C524A165D3EEAEDA68118F1EE7C17F3.
63.
Ge
H
,
Xu
K
,
Ghahramani
Z
.
Turing: a language for flexible probabilistic inference
.
In: International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9–11 April 2018
,
Playa Blanca, Lanzarote
,
Canary Islands, Spain
,
1682
90
(
2018
).
Available from:
http://proceedings.mlr.press/v84/ge18b.html.
64.
Faustino-Rocha
A
,
Oliveira
PA
,
Pinho-Oliveira
J
,
Teixeira-Guedes
C
,
Soares-Maia
R
,
da Costa
RG
, et al
.
Estimation of rat mammary tumor volume using caliper and ultrasonography measurements
.
Lab Anim
2013
;
42
:
217
24
.
65.
Del Monte
U
.
Does the cell number 109 still really fit one gram of tumor tissue?
Cell Cycle
2009
;
8
:
505
6
.
66.
Hoffman
MD
,
Gelman
A
.
The No-U-turn sampler: adaptively setting path lengths in hamiltonian monte carlo
.
J Mach Learn Res
2014
;
15
:
1593
623
.
67.
Blaom
A
,
Kiraly
F
,
Lienart
T
,
Simillides
Y
,
Arenas
D
,
Vollmer
S
.
MLJ: a julia package for composable machine learning
.
J Open Source Software
2020
;
5
:
2704
.
68.
Pedregosa
F
, et al
.
Scikit-learn: machine learning in Python
.
Journal of Machine Learning Research
2011
;
12
:
2825
30
.
69.
Hu
Z
,
Li
Z
,
Ma
Z
,
Curtis
C
.
Multi-cancer analysis of clonality and the timing of systemic spread in paired primary tumors and metastases
.
Nat Genet
2020
;
52
:
701
8
.
70.
Allen
LJS
,
van den Driessche
P
.
Relations between deterministic and stochastic thresholds for disease extinction in continuous- and discrete-time infectious disease models
.
Math Biosci
2013
;
243
:
99
108
.
71.
Kreger
J
,
Komarova
NL
,
Wodarz
D
.
A hybrid stochastic-deterministic approach to explore multiple infection and evolution in HIV
.
PLoS Comput Biol
2021
;
17
:
e1009713
.
72.
Tsimring
LS
.
Noise in Biology
.
Rep Prog Phys
2014
;
77
:
026601
.
73.
Movahedi
K
,
Guilliams
M
,
Van den Bossche
J
,
Van den Bergh
R
,
Gysemans
C
,
Beschin
A
, et al
.
Identification of discrete tumor-induced myeloid-derived suppressor cell subpopulations with distinct T cell-suppressive activity
.
Blood
2008
;
111
:
4233
44
.
74.
Luo
MC
,
Nikolopoulou
E
,
Gevertz
JL
.
From fitting the average to fitting the individual: a cautionary tale for mathematical modelers
.
Front Oncol
2022
;
12
:
793908
.
75.
Laghi
L
,
Bianchi
P
,
Miranda
E
,
Balladore
E
,
Pacetti
V
,
Grizzi
F
, et al
.
CD3+ cells at the invasive margin of deeply invading (pT3-T4) colorectal cancer and risk of post-surgical metastasis: a longitudinal study
.
Lancet Oncol
2009
;
10
:
877
84
.
76.
Galon
J
,
Costes
A
,
Sanchez-Cabo
F
,
Kirilovsky
A
,
Mlecnik
B
,
Lagorce-Pagès
C
, et al
.
Type, density, and location of immune cells within human colorectal tumors predict clinical outcome
.
Science
2006
;
313
:
1960
4
.
77.
Apodaca
MC
,
Wright
AE
,
Riggins
AM
,
Harris
WP
,
Yeung
RS
,
Yu
L
, et al
.
Characterization of a whole blood assay for quantifying myeloidderived suppressor cells
.
J Immunother Cancer
2019
;
7
:
230
.
78.
Peranzoni
E
,
Ingangi
V
,
Masetto
E
,
Pinton
L
,
Marigo
I
.
Myeloid cells as clinical biomarkers for immune checkpoint blockade
.
Front Immunol
2020
;
11
:
1590
.
79.
Bergenfelz
C
,
Larsson
A-M
,
von Stedingk
K
,
Gruvberger-Saal
S
,
Aaltonen
K
,
Jansson
S
, et al
.
Systemic monocytic-MDSCs are generated from monocytes and correlate with disease progression in breast cancer patients
.
PLoS One
2015
;
10
:
e0127028
.
80.
Benzekry
S
,
Tracz
A
,
Mastri
M
,
Corbelli
R
,
Barbolosi
D
,
Ebos
JML
.
Modeling spontaneous metastasis following surgery: an in vivo-in silico approach
.
Cancer Res
2016
;
76
:
535
47
.
81.
Sahoo
P
,
Yang
X
,
Abler
D
,
Maestrini
D
,
Adhikarla
V
,
Frankhouser
D
, et al
.
Mathematical deconvolution of CAR t-cell proliferation and exhaustion from real-time killing assay data
.
J R Soc, Interface
2020
;
17
:
20190734
.
82.
Rao
A
,
Barkley
D
,
França
GS
,
Yanai
I
.
Exploring tissue architecture using spatial transcriptomics
.
Nature
2021
;
596
:
211
20
.
83.
Andersson
A
,
Larsson
L
,
Stenbeck
L
,
Salmén
F
,
Ehinger
A
,
Wu
SZ
, et al
.
Spatial deconvolution of HER2-positive breast cancer delineates tumorassociated cell type interactions
.
Nat Commun
2021
;
12
:
6012
.
84.
Kuznetsov
VA
,
Makalkin
IA
,
Taylor
MA
,
Perelson
AS
.
Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis
.
Bull Math Biol
1994
;
56
:
295
321
.
85.
Diefenbach
A
,
Jensen
ER
,
Jamieson
AM
,
Raulet
DH
.
Rae1 and H60 ligands of the NKG2D receptor stimulate tumour immunity
.
Nature
2001
;
413
:
165
71
.
86.
Abedi-Valugerdi
M
,
Wolfsberger
J
,
Pillai
PR
,
Zheng
W
,
Sadeghi
B
,
Zhao
Y
, et al
.
Suppressive effects of low-dose 5-fluorouracil, busulfan or treosulfan on the expansion of circulatory neutrophils and myeloid derived immunosuppressor cells in tumor-bearing mice
.
Int Immunopharmacol
2016
;
40
:
41
49
.
87.
Pillay
J
,
den Braber
I
,
Vrisekoop
N
,
Kwast
LM
,
de Boer
RJ
,
Borghans
JAM
, et al
.
In vivo labeling with 2H2O reveals a human neutrophil lifespan of 5.4 days
.
Blood
2010
;
116
:
625
7
.
88.
Tak
T
,
Tesselaar
K
,
Pillay
J
,
Borghans
JAM
,
Koenderman
L
.
What's your age again? determination of human neutrophil half-lives revisited
.
J Leukocyte Biol
2013
;
94
:
595
601
.
89.
Yates
A
Callard
R
.
Cell death and the maintenance of immunological memory
.
Discrete & Continuous Dynamical Systems - B
2001
;
1
:
43
59
.
90.
Lanzavecchia
A
,
Sallusto
F
.
Dynamics of T lymphocyte responses: intermediates, effectors, and memory cells
.
Science
2000
;
290
:
92
97
.
91.
Huang
Y
,
Ma
C
,
Zhang
Q
,
Ye
J
,
Wang
F
,
Zhang
Y
, et al
.
CD4 + and CD8 + T cells have opposing roles in breast cancer progression and outcome
.
Oncotarget
2015
;
6
:
17462
78
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.