Abstract
Glioblastoma (GBM) is an aggressive primary brain cancer that currently has minimally effective treatments. Like other cancers, immunosuppression by the PD-L1-PD-1 immune checkpoint complex is a prominent axis by which glioma cells evade the immune system. Myeloid-derived suppressor cells (MDSCs), which are recruited to the glioma microenviroment, also contribute to the immunosuppressed GBM microenvironment by suppressing T cell functions. In this paper, we propose a GBM-specific tumor-immune ordinary differential equations model of glioma cells, T cells, and MDSCs to provide theoretical insights into the interactions between these cells. Equilibrium and stability analysis indicates that there are unique tumorous and tumor-free equilibria which are locally stable under certain conditions. Further, the tumor-free equilibrium is globally stable when T cell activation and the tumor kill rate by T cells overcome tumor growth, T cell inhibition by PD-L1-PD-1 and MDSCs, and the T cell death rate. Bifurcation analysis suggests that a treatment plan that includes surgical resection and therapeutics targeting immune suppression caused by the PD-L1-PD1 complex and MDSCs results in the system tending to the tumor-free equilibrium. Using a set of preclinical experimental data, we implement the approximate Bayesian computation (ABC) rejection method to construct probability density distributions that estimate model parameters. These distributions inform an appropriate search curve for global sensitivity analysis using the extended fourier amplitude sensitivity test. Sensitivity results combined with the ABC method suggest that parameter interaction is occurring between the drivers of tumor burden, which are the tumor growth rate and carrying capacity as well as the tumor kill rate by T cells, and the two modeled forms of immunosuppression, PD-L1-PD-1 immune checkpoint and MDSC suppression of T cells. Thus, treatment with an immune checkpoint inhibitor in combination with a therapeutic targeting the inhibitory mechanisms of MDSCs should be explored.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
Glioblastoma (GBM) is the most common and aggressive type of primary brain cancer, claiming tens of thousands of lives each year. Within 5 years of diagnosis, less than 10% of GBM patients survive, and most succumb to the tumor within 15 months of diagnosis—even with treatment (Ostrum et al. 2015, Fernandes et al. 2017). The current standard of care includes surgical resection followed by radiotherapy and chemotherapy with temozolomide (TMZ) (Stupp et al. 2005, Fernandes et al. 2017), but the minimally effective results demonstrate a need to look for new treatment strategies. The highly complex and immune suppressed tumor microenvironment of GBM makes treatment difficult, thus immunotherapy, which has shown to be effective in some cancer types, holds promise (Brown et al. 2018).
In the last decade, an increasing number of immunotherapies have been developed for GBM (Yu and Quail 2021, Bausart et al. 2022, Bryukhovetskiy 2022). Anti-PD-1 is a standard immunotherapy that has proven to be beneficial in treating other cancers, such as melanoma (Postow et al. 2015), Hodgkin’s lymphoma (Ansell et al. 2015), colon (Duraiswamy et al. 2013), cervical (Chung et al. 2019), and non-small cell lung cancer (NSCLC) (Brahmer et al. 2015), although it has failed as a monotherapy in GBM phase III clinical trials (Preusser et al. 2015, Reardon et al. 2017, Lim et al. 2018). Cancer can evade the immune system by expressing the marker PD-L1, which downregulates the cytotoxic response of activated T cells by binding to their PD-1 receptor. Treatment with anti-PD-1 unmasks the tumor by binding to PD-1, thereby facilitating tumor recognition by T cells and enhancing the immune response. The failure of anti-PD-1 monotherapy in GBM might potentially be explained by additional mechanisms of immune suppression, such as the infiltration of myeloid-derived suppressor cells (MDSCs), which are cells that suppress T cell activity.
Gliomas mediate recruitment of CCR2\(^+\) MDSCs by releasing chemokines such as CCL2 and CCL7, which are cognate ligands of the CCR2 receptor (Takacs et al. 2021, Takacs et al. 2022). Higher expression of CCL2 and CCL7, as well as CCR2, within the tumor microenvironment all correlate with a worse prognosis for gliomas (Korbecki et al. 2020, Chang et al. 2016). Once the MDSCs enter the tumor microenvironment, they suppress T cells through a variety of mechanisms including, but not limited to, inducing apoptosis in activated T cells (Saio et al. 2001), producing enzymes (such as arginase) which metabolize amino acids needed for T cell proliferation (Condamine and Gabrilovich 2011), and secreting other factors which modulate the immune response (NO, TGF-\(\beta \), TNF-\(\alpha \), H\(_2\)O\(_2\)) (Monu and Frey 2012, Markowitz et al. 2017). These mechanisms, along with MDSC promotion of angiogenesis (Vetsika et al. 2019), lead to rapid tumor progression. Treating with a CCR2 antagonist aids the immune response by reducing MDSC recruitment and thus sequestering the MDSC population to the bone marrow. Flores-Toro et al. (2020) reported improved survival of two murine glioma models (KR158 and 005 GSC) with combination treatment of a CCR2 antagonist (CCX872) and anti-PD-1.
Although MDSCs play a significant role in glioma progression, there are few published mathematical models of tumor–immune dynamics that incorporate these cells. Allahverdy et al. (2019) presents a discrete agent-based model of MDSCs in the context of a general tumor and simulates treatment with the chemotherapeutic agent 5-fluorouracil (5-FU). Similar to Allahverdy et al. (2019), the ODE model of Shariatpanahi et al. (2018) for a general tumor also simulates treatment with 5-FU. Further, they include treatment with L-arginine since the model considers production of arginase I as the primary mechanism by which MDSCs suppress T cells. Liao et al. (2014) highlights a mechanism of MDSC recruitment by developing a PDE model which focuses on the role of interleukin-35 (IL-35) in MDSC recruitment and tumor growth. Lai et al. (2018) used a PDE model of breast cancer to determine the influence of the immune checkpoint inhibitor, anti-CTLA-4, on M2 macrophages within the tumor site. Although they did not specifically model MDSCs, their work implied that these conclusions apply to MDSCs because a subset of monocytic MDSCs differentiate into M2 macrophages. Kreger et al. (2023) developed a stochastic delay differential equations model of a metastasizing breast cancer to evaluate MDSC influence on tumor progression. Although the model is fitted to patient response data to immune checkpoint inhibitors, the model itself does not include a mechanism to represent an immune checkpoint. Each of these models include MDSCs and provide a foundation to incorporate the PD-L1-PD-1 immune checkpoint. Furthermore, we aim to develop a model specifically for GBM instead of a general tumor.
Mathematical models of tumor–immune dynamics range from modeling general tumors (Lai and Friedman 2017, Nikolopoulou et al. 2018, Eftimie et al. 2011, Mahlbacher et al. 2019, Shi et al. 2021, Radunskaya et al. 2018, Khyat and Jang 2022) to specific cancers including GBM, non-small lung cancer (NSCLC), melanoma, and breast cancer (Storey et al. 2020, Yu and Jang 2019, Butner et al. 2021, Lai et al. 2018, Banerjee et al. 2015, Özköse et al. 2022, Mirzaei et al. 2021, Khajanchi and Banerjee 2017, Khajanchi 2021, Jafarnejad et al. 2019, Perlstein et al. 2019, Bitsouni and Tsilidis 2022). Similar to Lai et al. (2018), Yu and Jang (2019) modelled the CTLA-4 immune checkpoint and included anti-CTLA-4 therapy, but restricted their focus to a four equation ODE model of tumor cells, CD4\(^+\) T cells, IFN-\(\gamma \), and CTLA-4. Two years earlier, Lai and Friedman (2017) modelled a different immune checkpoint, the PD-L1-PD-1 complex, by developing a PDE model of a general tumor which included additional cells and molecules like dendritic cells, CD4\(^+\)/CD8\(^+\) T cells, IL-2, and IL-12. Their model was used to simulate treatment with anti-PD-1 and a cancer vaccine, GVAX. This work was later reduced to an ODE model by Nikolopoulou et al. (2018) which focused on tumor and T cells with intermittent and continuous anti-PD-1 treatment. Shi et al. (2021) and Nikolopoulou et al. (2021) extended this work by establishing the global dynamics of the model along with a generalized version. Storey et al. (2020) applied these models to GBM by estimating parameters from in vivo murine experimental data. Their ODE model divided the immune compartment into the categories of innate and adaptive immune cells and modeled treatment with anti-PD-1 and an oncolytic viral therapy. Radunskaya et al. (2018) deviated from modeling the typical immune checkpoint inhibitor, anti-PD-1, by considering treatment with anti-PD-L1 in their model of general tumor-immune dynamics, which compartmentalized interactions to be within the spleen, blood, and tumor site. While each of these previous papers focused on a specific immune checkpoint, Butner et al. (2021) developed an ODE model which was not necessarily specific to the type of immune checkpoint or the type of tumor. Using this model along with patient data, Butner et al. (2021) predicted the long-term tumor burden for a variety of patients being treated with different immune checkpoint inhibitors.
In this paper, we develop a GBM-specific ODE model of cancer and T cell interactions by incorporating the PD-L1-PD-1 immune checkpoint along with tumor recruitment of MDSCs. Parameter values are estimated from the literature and via comparison with experimental in vivo murine data.
Our paper unfolds as follows: in Sect. 2, we describe the GBM-immune dynamics model. In Sect. 3, we perform equilibrium and solution analysis and emphasize these findings through bifurcation analysis. Section 4 analyzes the parameter space by implementing the Approximate Bayesian Computation (ABC) method used in conjunction with the extended Fourier Amplitude Sensitivity Test (eFAST). Numerical simulations alongside experimental data are also provided. In Sect. 5, we conclude with a discussion of our results and future directions.
2 Mathematical model
The glioblastoma (GBM)–immune model focuses on the dynamics within the glioma microenvironment and is based on the previous work of Nikolopoulou et��al. (2018), Storey et al. (2020), and Shariatpanahi et al. (2018). A diagram of the biological interactions implemented is illustrated in Fig. 1, and Table 1 lists the biological meaning of the model’s parameters with their respective units and representative range of values.
Letting C be the number of tumor cells, T the number of activated T cells (all cells that have the CD3\(^+\) marker, which includes both CD4\(^+\) and CD8\(^+\) T cell subsets), and M be the number of myeloid-derived suppressor cells (MDSCs), we arrive at the following system of equations:
In Eq. (1a), we assume the tumor cells exhibit logistic growth with a carrying capacity of \(C_\textrm{max}\). T cells, which are the main drivers of tumor cell death, kill glioma cells at a rate of \(\eta \), as in Nikolopoulou et al. (2018).
Equation (1b) is of the same form as in Lai and Friedman (2017), Nikolopoulou et al. (2018), and Storey et al. (2020). The first term represents T cells activation at a constant rate \(a_T\), potentially due to the cytokine IL-12 which is explicitly represented in the previous models, and the recruitment (or stimulation) of T cells to the tumor site by the presence of a glioma at a rate \(s_T\). Once in the glioma microenvironment, T cells are exhausted by the formation of the PD-L1-PD-1 complex, which binds PD-L1 present on glioma cells and T cells with PD-1 on T cells. Assuming \((T + \epsilon _C C)\) characterizes the level of free PD-L1 in the tumor microenvironment, the PD-L1-PD-1 complex is represented by \(T(T + \epsilon _C C)\) since only T cells express PD-1. As T cells and tumor cells increase, more PD-L1-PD-1 complexes form, thus decreasing the overall fraction to represent T cell inhibition. The parameter \(\rho \) results from combining multiple parameters in Lai and Friedman (2017), Nikolopoulou et al. (2018), and Storey et al. (2020) due to their structural non-identifiability. The remaining terms of (1b) represent T cells suppressed/deactivated by MDSCs (Gabrilovich and Nagaraj 2009) by a rate r and T cells dying at a rate \(d_T\).
In Eq. (1c), the first term represents glioma recruitment of MDSCs to the microenvironment by secreting the chemokines CCL2 and CCL7, which are ligands for the CCR2 receptor expressed by MDSCs. This results in a chemotactic effect, drawing MDSCs from the bone marrow to the glioma microenvironment, where here we are specifically focusing on monocytic M-MDSCs as opposed to granulocytic PMN-MDSCs. Since these chemokines are produced by glioma cells, this recruitment increases as glioma cells increase. We assume that the presence of a tumor also results in the expansion of splenic MDSCs (Liu et al. 2007, Fig. 2), using the same fractional term as the model of Shariatpanahi et al. (2018) in the second term of (1c), thus indirectly causing more MDSCs to accumulate in the tumor. The last term represents MDSCs dying at a rate \(d_M\).
Our model is distinct from previous models in that it combines both the PD-L1-PD-1 immune checkpoint, such as incorporated in Lai and Friedman (2017), Nikolopoulou et al. (2018), and Storey et al. (2020), and MDSCs, such as in Shariatpanahi et al. (2018), and it additionally was developed specifically for GBM. Though Storey et al. (2020) also focused on GBM, while the other studies were for general tumors or tumors different than GBM, their study included oncolytic viral therapy treatment. Shariatpanahi et al. (2018) differs from our model in that MDSC dynamics are modeled in the spleen instead of the tumor site and they incorporated chemotherapy treatment. Here, we integrate an immune checkpoint mechanism and MDSCs to study each of these immunotherapy targets in GBM.
3 Equilibrium and solution analysis
We describe the equilibria and stability of system (1) in this section to verify that our model is biologically reasonable and to determine suitable target parameters for treatment. Since the functions on the right hand side of the equations in system (1) are continuously differentiable, solutions exist and are unique by standard ordinary differential equations theory.
We begin by confirming that our system will neither veer into negative cell counts nor increase to infinity over time but will be realistically bounded. First, we show that solutions with positive initial conditions remain positive for all time.
Theorem 1
(Positivity) Solutions of (1) that start positive remain positive.
Proof
Assume \(C(t_0)\) > 0, \(T(t_0)\) > 0, and \(M(t_0)\) > 0 for some initial time \(t_0\). Without loss of generality, assume \(t_0 = 0\).
In order to have T(t) \(\le \) 0 for some \(t>0\), it is necessary to have dT/dt \(\le \) 0 for \(T=0\). However, when \(T=0\), (1b) becomes \(dT/dt =a_T >0\). Thus, T(t) > 0 for all \(t>0\).
Assuming that C and T exist on some interval [0, t) for \(t>0\), when we integrate (1a), we arrive at
which is always positive since C(0) > 0.
Lastly, when \(M=0\), (1c) becomes
Since we proved that C(0) > 0 implies C(t) > 0 for all \(t>0\) and because all the parameters are positive, we conclude that for \(M=0\), \(dM/dt > 0\). Thus, M(t) > 0 for all \(t>0\). \(\square \)
Next, we show that solutions of system (1) with positive initial conditions are bounded.
Theorem 2
(Boundedness) Solutions of (1) that start positive are bounded.
Proof
Assume \(C(t_0)\) > 0, \(T(t_0)\) > 0, and \(M(t_0)\) > 0 for some initial time \(t_0\). Without loss of generality, assume \(t_0 = 0\).
Since all parameters are positive and C and T are positive by Theorem 1, by comparison we can bound (1a) by
Since logistic growth is bounded, C(t) is bounded as well. Specifically,
for all \(t\ge 0\).
For the T cells, note that (1b) can be bounded as
g(T) achieves an absolute maximum at
Hence, (6) can be further bounded as
Using a standard comparison argument,
Thus, for all \(t\ge 0\),
For the MDSCs, since C is bounded (5), (1c) can be bounded as
By a standard comparison argument,
Therefore, for all \(t\ge 0\),
Hence, system (1) is bounded. \(\square \)
The system (1) admits two categories of fixed points: tumor-free equilibria of the form \((C,T,M)=(0,T_0^*,M_0^*)\) and tumorous equilibria of the form \((C,T,M)=(C^*,T^*,M^*)\). We first show that there exists one unique tumor-free equilibrium.
Theorem 3
(Uniqueness of tumor-free equilibrium) The system (1) has a unique tumor-free equilibrium \((C,T,M) = (0,T_0^*,0)\), where \(T_0^*\) is increasing with respect to \(a_T\) and decreasing with respect to \(d_T\) and \(\rho \).
Proof
Assume that there is no tumor present, i.e., \(C =0\). Setting \(dM/dt=0\) in (1c) implies that \(M_0^* = 0\).
When \(dT/dt = 0\), (1b) implies that \(d_T T_0^* (1+\rho {T_0^*}^2) =a_T\). Let
According to Descartes’ rule of signs, the number of sign changes in front of the coefficients of a polynomial corresponds to the number of positive zeros. Since there is exactly one sign change because all the parameters are positive, there is exactly one positive zero, \(T_0^*\). Thus, a biologically relevant tumor-free equilibrium \((0,T_0^*,0)\) exists and is unique.
Note that as \(a_T\) increases, the graph of f(T) shifts down, resulting in a larger \(T_0^*\). As \(d_T\) and \(\rho \) increase, \(f'(T)=2d_T \rho T^2 +d_T\) increases for \(T>0\), which results in a narrowed curve and thus a smaller \(T_0^*\).
Lastly, when we solve for \(T_0^*\) using the cubic formula, we obtain
\(\square \)
Next, we determine conditions for the stability of the tumor-free equilibrium \((0,T_0^*,0)\).
Theorem 4
(Stability of tumor-free equilibrium) When \(\lambda _C < \eta T_0^*\), the tumor-free equilibrium, \((0,T_0^*,0)\), is locally asymptotically stable. However, when \(\lambda _C > \eta T_0^*\), \((0,T_0^*,0)\) is a saddle point.
Proof
The Jacobian evaluated at the tumor-free equilibrium is
and the eigenvalues are
Therefore, when \(\lambda _C < \eta T_0^*\), the tumor-free equilibrium \((0,T_0^*,0)\) is locally asymptotically stable. However, when \(\lambda _C > \eta T_0^*\), \(\lambda _1>0\), so the tumor-free equilibrium is a saddle point. \(\square \)
Biologically, the condition \(\lambda _C<\eta T_0^*\) guaranteeing local stability of the tumor-free equilibrium indicates that T cells kill cancer cells faster than the cancer cells can multiply during the onset of tumor initiation. As for treatment, the worst patient scenario would be a saddle tumor-free equilibrium, as this would suggest that the patient will most likely relapse regardless of their proximity to the disease-free condition. Thus, by varying \(\lambda _C\), \(\eta \), or \(T_0^*\) (15) (which depends on \(a_T\), \(d_T\), and \(\rho \)), a patient’s tumor-free equilibrium could shift from a saddle to a locally stable point. These parameters would likely be affected in an immune-compromised individual, and in particular, treatment with an immune checkpoint inhibitor to decrease \(\rho \) could be an effective strategy.
In Appendix A, we establish conditions for the uniqueness (Theorem 6) and local stability (Theorem 7) of the tumorous equilibrium. We conclude our discussion here by determining conditions for the global stability of the tumor-free equilibrium \((0,T_0^*,0)\).
Theorem 5
(Global stability of tumor-free equilibrium) The system (1) has a globally stable tumor-free equilibrium when \(\lambda _C < \eta \beta \), where \(\beta := a_T /(1+\rho \widehat{T}(\widehat{T}+ \epsilon _C \widehat{C}))(r\widehat{M} +d_T)\) and \(\widehat{C}, \ \widehat{T}, \ \text {and} \ \widehat{M}\) are the upper bounds for glioma cells (5), T cells (10), and MDSCs (13), respectively, as determined in Theorem 2.
Proof
From Theorems 1–2, we have that C(t), T(t), and M(t) are bounded above with upper bounds denoted by \(\widehat{C}\), \(\widehat{T}\), and \(\widehat{M} >0\), respectively, and below by 0. By these bounds,
Thus, using a standard comparison argument,
If \(\lambda _C < \eta \beta \), then for all \(s>0\) sufficiently large,
which, taking the limit of (2), implies that
Bounding T(t) in the tumor cell equation (1a), we have
which, taking the limit as \(t\rightarrow \infty \), implies that \(dC/dt \rightarrow 0\). Now, using (21),
where \(\epsilon \rightarrow 0\) as \(t \rightarrow \infty \). A standard comparison argument reveals that
for any \(\epsilon >0\). Thus, since M(t) is nonnegative for all \(t\ge 0\), \(M(t) \rightarrow 0\) as \(t \rightarrow \infty \), and hence, \(dM/dt \rightarrow 0\) as \(t\rightarrow \infty \).
Since T is bounded above and below and since \(C(t) \rightarrow 0\) and \(M(t) \rightarrow 0\) as \(t\rightarrow \infty \), the limiting equation of (1b) is \(dT/dt=G(T)\), where
and we note that \(G(T_0^*)=0\). Since
where \(\epsilon _t \rightarrow 0\) as \(t \rightarrow \infty \), using a comparison argument,
Thus, \(T(t) \rightarrow T_0^*\) as \(t\rightarrow \infty \), and \(dT/dt \rightarrow G(T_0^*) = 0\) as \(t \rightarrow \infty \).
This allows us to conclude that the tumor-free equilibrium point, \((0,T_0^*,0)\), is globally stable in \(\mathbb {R}_0^+ \times \mathbb {R}_0^+\) under the condition that \(\lambda _C < \eta \beta \). \(\square \)
The best patient scenario is one in which the tumor will tend toward the tumor-free equilibrium regardless of size. While satisfying the global stability condition in Theorem 5 is unlikely in reality, as we shall see in Sect. 4.3, the condition provides additional targets for treatment, such as \(\epsilon _C\), r, and \(C_\textrm{max}\). These parameters could be targeted through treatment with anti-PD-1/PD-L1, L-arginine, and surgical resection, respectively. Through parameter analysis in Sect. 4, we will narrow our search for therapeutic targets.
3.1 Bifurcation analysis
We validate the tumor-free equilibrium global stability conditions (Theorem 5) with a bifurcation analysis. Figure 2 indicates that the tumor growth rate (\(\lambda _C\)), T cell kill rate of tumor cells (\(\eta \)), and T cell inhibition rate by PD-L1-PD-1 (\(\rho \)) and MDSCs (r) are reasonable targets for treatment, since shifting these parameters can result in the system tending toward the tumor-free equilibrium. Varying the remaining parameters cannot considerably relieve the tumor burden (Fig. S1 (Online Resource 1, Sect. 5)) and thus are unlikely to be advantageous treatment targets.
Regarding MDSC infiltration and suppression, Fig. 2 suggests that targeting MDSC recruitment, death, and their inhibition of T cells can increase the cytotoxic immune response; however, only by directly targeting immune suppression by MDSCs, r, can tumor cells be reduced to zero. Thus, therapeutics which prevent T cell inhibition by MDSCs should be used to relieve the tumor burden. These should be used in combination with immunotherapies targeting T cell inhibition by PD-L1-PD-1 (\(\rho \)), such as anti-PD-1 or anti-PD-L1, since varying \(\rho \) also minimized the tumor size (Fig. 2).
4 Parameter analysis
In this section, we analyze the parameter space and its effect on system (1) by implementing the Approximate Bayesian Computation (ABC) rejection method (Sect. 4.2) to compare with experimental data (Sect. 4.1), plotting numerical simulations (Sect. 4.3), and performing global sensitivity analysis using the extended Fourier Analysis Sensitivity Test (eFAST) (Sect. 4.4). Section 4.5 contains an examination of the results.
4.1 Experimental data
The model of GBM–immune dynamics (1) developed in this study is applied to a set of experimental data from a high-grade murine glioma cell line, KR158.
Mice were anesthetized and surgically implanted with 35,000 glioma cells through an incision in the skull. At several time points (7, 13, 20, 24, 27, and 34 days after glioma cell implantation), 1–4 mice were euthanized followed by resection of the brain. 2D images of the glioma tissue were used to determine tumor, T cell, and MDSC cell counts (Fig. 3 and Online Resource 2).
For a more detailed description of the experimental setup, image analysis pipeline, and data conversion from a 2D section to a 3D spherical tumor, see Online Resource 1 (Sects. 1 and 2).
4.2 Approximate Bayesian computation
Deterministic approaches, like the gradient descent method, are commonly used to estimate a single set of parameter values which optimally represent the given data (Allmaras et al. 2013). This optimization is computed by minimizing error between the numerical simulation and data points. However, depending on the initial guesses of these parameters, this method can converge to different local error minima, making it difficult to uncover a set of parameter estimates which is globally optimal. Additionally, these approaches do not quantify uncertainty in the parameter estimation and generally may not give reliable fits for small data sets. To remedy these issues, we utilize the Approximate Bayesian Computation (ABC) rejection method (Sunnåker et al. 2013, Liepe et al. 2014). This algorithm produces probability distributions of the parameter values, which can indicate a range of near-optimal values for each parameter. These distributions provide more detailed information should additional data be made available for a more precise data fitting. In recent years, the use of the ABC rejection method in the study of complex biological systems has been increasing (Browning et al. 2017, Stepien et al. 2019, Xiao et al. 2021, da Costa et al. 2018).
The first step of the ABC rejection method is to specify a prior distribution for each of the parameters. Consequently, we generate 1 million parameter sets, i.e., 1 million 13-dimensional vectors, where each vector forms a parameter set
The values for each parameter were chosen uniformly from the ranges shown as the bounds of the horizontal axes in Fig. 4.
Next, we calculate the relative error of a single parameter set,
where \(d_i\) is the experimental data point and \(x_i\) is the simulated value at the \(i^\textrm{th}\) time point. We set \(d_i\) to be the average cell count for the \(i^\textrm{th}\) time point, where the time points were days 7, 13, 20, 24, 27, and 34. The relative errors for tumor cells (\(E_C\)), T cells (\(E_T\)), and MDSCs (\(E_M\)) are calculated individually and the total error is \(E_\textrm{total} = E_C+E_T+E_M\).
Lastly, we specify an error threshold, R: if the error of a parameter set is less than the threshold, the set is accepted (i.e., the simulation is sufficiently close to the experimental data), but if it is greater than the threshold, the parameter set is rejected. Since we observed that parameter sets that matched the tumor cell and MDSC data well often did not match the T cell data, we set error thresholds separately for each cell type. In particular, we set \(R_C=0.75\) for the tumor cell error threshold, \(R_T=0.72\) for the T cells, and \(R_M=0.78\) for the MDSCs. For a parameter set to be accepted, we required \(E_C\le R_C\), \(E_T\le R_T\), and \(E_M\le R_M\).
Of the 1 million parameter sets sampled, approximately 44,000 sets were accepted. Visualizations of the resulting posterior distributions are shown in Fig. 4. Along the diagonal of Fig. 4a are smoothed 1D histograms for parameters exhibiting a right-skewed distribution, while Fig. 4b shows the remaining parameters. Below the diagonal of both of these figures are 2D projections of pairs of parameters. Red shaded areas correspond with higher frequency of the parameter in the posterior distribution.
Table 2 lists the summary statistics of the posterior distributions as well as the best fitting probability distribution, which was calculated by minimizing the 1-Wasserstein metric (earth mover’s distance).
To determine whether our results were robust, we extracted different percentages of the accepted parameters and compared parameter distributions. Distributions remained largely unchanged when we analyzed any percentage of our accepted parameters.
Since the posterior distributions for parameters in Fig. 4b have less distinct peaks (cf. Fig. 4a), this suggests that, given our data, we were unable to attain additional information on the actual probability distributions for these parameters beyond upper and lower bounds found in the literature.
4.3 Numerical simulations
Since mice were orthotopically implanted with 35,000 KR158 glioma cells in the experiments (Sect. 4.1), myeloid-derived suppressor cells (MDSCs) are absent from the brain before tumor introduction, and the number of activated T cells are initially negligible, the initial conditions for numerical simulations were set to be \(C(0) = 35,000\), \(T(0) = 0\), and \(M(0) = 0\) cells. We assume that any loss of glioma cells due to implantation is negligible.
We ran 10,000 simulations using parameter sets sampled from the ABC method posterior distributions listed in Table 2, and illustrate the mean cell counts and cell counts within standard deviation \(\sigma /4\), \(\sigma /2\), \(3\sigma /4\), and \(\sigma \) calculated at each hour in Fig. 5a. We additionally plot the experimental data points for glioma cells, T cells, and MDSCs and observe that the simulations capture nearly all data points within one standard deviation from the mean.
The ABC rejection method histograms illustrated in Fig. 4 were produced by calculating the relative error (28), where \(d_i\) are the cell counts averaged at each time point. However, we also used the minimum and maximum cell counts for \(d_i\) to evaluate the relative error in two separate instances of the ABC method. The parameter sets that gave rise to minimal error from these three different trials are listed in Table S1 (Online Resource 1, Sect. 4) and numerical simulations are shown in Fig. 5b. A comparison of their total error values could suggest that our parameter ranges are better suited for larger gliomas (error: 0.937) rather than smaller gliomas (error: 1.78), while they are best suited for the average glioma size (error: 0.785). However, these differences in error could also be explained by the variability in the data.
We find that the parameter sets from these three trials (Table S1 in the Online Resource 1, Sect. 4) each satisfy the conditions for a saddle tumor-free equilibrium (Theorem 4), thus corresponding to the worst treatment scenario. We cannot guarantee the existence of a tumorous equilibrium, \((C^*, T^*, M^*)\), for the minimum and average data parameter sets since they fail to fulfill \(C_\textrm{max} \epsilon _C \eta / \lambda _C <1\) in Theorem 6 given in Appendix A. However, the maximum data parameter set exhibits a tumorous equilibrium when \(T^* < 5.92 \times 10^8 = \lambda _C / \eta \). Further, if this equilibrium exists, it is guaranteed to be unique and it is locally asymptotically stable when \(T^* \in (11, \ 2.96 \times 10^8)\). All T cell data points are contained within this range, so it is likely that large tumors have a unique, locally asymptotically stable tumorous equilibrium, making these tumors difficult to treat as the system tends towards a tumorous state.
4.4 Sensitivity analysis
To determine the influence of parameters on tumor progression, we conducted a sensitivity analysis using the extended Fourier Analysis Sensitivity Test (eFAST) method outlined in Saltelli et al. (1999) by utilizing the MATLAB code developed by the Kirshner Lab at the University of Michigan (Kirschner and Panetta 1998, Kirschner 2008). eFAST is a global sensitivity analysis method which neither requires model linearity nor a monotonic relationship between the output and the parameters (Saltelli et al. 2008). Similar to the Sobol’ method, it is a variance-based method, but it uses a monodimensional Fourier expansion of the model to evaluate variance. Because of this monodimensional transformation, each parameter i can be sampled within the unit hypercube by varying \(s \in (-\pi ,\pi )\) along the space-filling search curve, \(x_i\), defined by
where \(\omega _i\) is the angular frequency assigned to parameter i, \(\varphi _i\) is a random phase-shift, and \(G_i\) is a transformation function.
The ABC posterior distributions of each parameter (Fig. 4 and Table 2) are used to determine the transformation function, \(G_i\), hence the search curve, \(x_i\). For parameters with right-skewed posterior distributions (Fig. 4a: \(\lambda _C\), \(C_\textrm{max}\), \(\eta \), \(\rho \), \(\epsilon _C\), r), we used the search curve proposed by Cukier et al. (1973),
where \(n_i\) is the nominal value of the \(i\textrm{th}\) parameter, set to be the mode of the ABC posterior distribution, and \(v_i\) accounts for uncertainty by altering the width of the parameter range. \(v_i\) was set to be such that parameters are sampled within the interval [0, 1] by setting \(v_i = \ln (1/n_i)\).
For the remaining parameters (Fig. 4b: \(a_T\), \(s_T\), \(d_T\), \(s_M\), \(\alpha \), q, \(d_M\)), we used the search curve for a uniform distribution that was determined by Saltelli et al. (1999),
Using these search curves, we generated 2,049 parameter sets for each 5 resamplings, which resulted in a total sample size of 10,245 parameter sets. Details on how we chose the sample size and number of resamplings are in Online Resource 1 (Sect. 6.1). After parameters were sampled from the unit hypercube, the samples for each parameter were translated to the ranges shown as the bounds of the horizontal axes in Fig. 4.
The main effect, or first-order index, \(S_i\), describes the effect that a single parameter i has on the variance of the system independently of other parameters. In contrast, the total effect index, \(S_{Ti}\), is a sum of \(S_i\) and the higher-order interactions of parameter i with other parameters. In Fig. 6, we illustrate the main effect \(S_i\) and the total effect index \(S_{Ti}\) for each parameter calculated after a simulation time of 40 days (the approximate time at which an untreated mouse requires euthanasia). The larger the difference between \(S_{Ti}\) and \(S_i\), the more interaction there is between i and other parameters in synergistically influencing the variance of the output.
Figure 6 indicates that T cells are sensitive to perturbations in several parameters, including the kill rate of tumor cells by T cells (\(\eta \)) as well as the inhibition of T cells by PD-L1-PD-1 (\(\rho \)) and MDSCs (r). Tumor cells express sensitivity to the tumor growth rate (\(\lambda _C\)), tumor carrying capacity (\(C_\textrm{max}\)), and the kill rate by T cells (\(\eta \)). These results are similar for the MDSC population, however, the death rate (\(d_M\)) holds a significantly larger influence on variations. Also, while MDSCs and tumor cells display similar total order indices for \(\eta \) and \(\lambda _C\), the first order indices of these parameters are larger for tumor cells than MDSCs. This is most likely since these parameters directly influence the tumor cell population but indirectly affect the MDSC population.
Additional sensitivity analyses were run with simulation ending times of 5, 10, 20, 60, 80, and 100 days. Results for days 60, 80, and 100 were qualitatively similar to Fig. 6, and results for days 5, 10, and 20 are shown in Fig. S2 in Online Resource 1 (Sect. 6.2). Comparing the outputs, we find that each cell population becomes increasingly sensitive to \(C_\textrm{max}\) as time progresses, as limited space and nutrients become more significant for larger tumors. MDSCs are initially very sensitive to the MDSC recruitment rate (\(s_M\)), but this decreased with time and, by day 40, its influence is surpassed by \(d_M\). This result suggests that a CCR2 antagonist should be used as soon as possible, and if the tumor is in its later stages, therapeutics that kill MDSCs at the tumor site should be considered instead. T cells initially were the most sensitive to \(\rho \) but this switched to \(\eta \) by day 10. Further, T cell sensitivity to r continues to increase until it surpasses \(\rho \) by day 20, where it remains the second most influential parameter. Therefore, while treatment with an immune checkpoint inhibitor, such as anti-PD-1 or anti-PD-L1, should be beneficial at any time, results indicate that it is best to treat sooner if possible and to incorporate treatment strategies targeting T cell inhibition by MDSCs later in treatment.
Since the search curve (30) of Cukier et al. (1973) produces a stronger skew than our ABC distributions suggest (Fig. 4), we additionally sampled all parameters uniformly using (31) (Fig. S3 in Online Resource 1, Sect. 6.2). We assume that the actual sensitivity indices are within the range produced by the results in Figs. 6 and S3. Figure S3 shows slightly more interaction occurring between parameters (as seen by greater differences between \(S_i\) and \(S_{Ti}\)), thus, by utilizing the ABC posterior distributions to inform the search curve, Fig. 6 more precisely determined which parameters should be targeted with treatment.
4.5 Results
In the existing literature, we could find no estimates for the value of \(s_M\). From the ABC posterior distribution, we expect \(s_M\) to fall within the range of 0.0214–0.0770 day\(^{-1}\) (Table 2). Separately, we approximated \(s_M\) using data on CCL2 expression in KR158 glioma cells as well as the migratory effect of CCL2 on MDSCs in vitro (Online Resource 1, Sect. 3.1) and found that KR158 gliomas recruit MDSCs at a rate of 0.0123–0.0144 MDSCs per day per tumor cell. Since this data only included the effect of CCL2 and not CCL7 on MDSC migration, we expect this range to be an underestimate for \(s_M\). However, as the mode of the \(s_M\) distribution, 0.0249, is close to double the experimentally calculated value, we expect that CCL7 will have a similar effect on MDSC migration as CCL2, assuming that the interaction between CCL2 and CCL7 causes an additive effect in recruitment.
In all three populations, the tumor growth rate (\(\lambda _C\)), tumor carrying capacity (\(C_\textrm{max}\)), and the tumor kill rate by T cells (\(\eta \)) are sensitive parameters which exhibit interactions with other parameters, as seen by large differences between the first and total order indices. Although the sensitivity analysis indicates this interaction (Fig. 6), it does not identify which parameters are interacting with each other. By combining the sensitivity analysis results with the ABC histograms in Fig. 4a, we can conjecture potential relationships.
There is a considerable inverse relationship between \(\lambda _C\) and the PD-L1-PD-1 parameters (\(\rho \) and \(\epsilon _C\)), \(C_\textrm{max}\) and the inhibition by PD-L1-PD-1 (\(\rho \) and \(\epsilon _C\)) and MDSCs (r), as well as \(\eta \) and r and slightly between \(\eta \) and \(\epsilon _C\). Although these inverse relationships could potentially be due to these parameters each exhibiting right-skewed distributions, the inverse relationships between \(\lambda _C\), \(C_\textrm{max}\), and \(\eta \) themselves are negligible. Therefore, we hypothesize that, although the tumor population is less sensitive to the immune suppression parameters \(\rho \), \(\epsilon _C\), and r, targeting these would affect the tumor’s sensitivity to \(\lambda _C\), \(C_\textrm{max}\), and \(\eta \). At the very least, according to the T cell sensitivity in Fig. 6, targeting these parameters would increase the immune response, thus decreasing the tumor population.
Fig. 4a also displays an inverse relationship between the tumor upregulation of PD-L1 (\(\epsilon _C\)) and the inhibition of T cells by the PD-L1-PD-1 complex (\(\rho \)) and by MDSCs (r). The relationship between \(\rho \) and \(\epsilon _C\) suggests that as tumor cells express more PD-L1, the PD-L1-PD-1 complex need not be as effective at inhibition in order to produce the same outcome. The second relationship suggests that r, although it is a less sensitive parameter, helps to increase the sensitivity of the system to \(\epsilon _C\). Thus, perturbations in the two forms of immunosuppression, the PD-L1-PD-1 complex and MDSCs, cause greater variance on the tumor together rather than alone.
5 Discussion and future directions
Current treatment paradigms for the highly aggressive brain tumor glioblastoma (GBM) present limited efficacy. Barriers, such as the highly immune suppressed tumor microenvironment, make gliomas difficult to treat. To address this barrier, we mathematically model GBM-immune dynamics by considering the role of myeloid-derived suppressor cells (MDSCs) in aiding GBM progression through T cell suppression. Gliomas recruit CCR2\(^+\) MDSCs to the tumor microenvironment by expressing the chemokines CCL2 and CCL7, which are ligands of the CCR2 receptor. Once at the tumor site, MDSCs suppress T cells, which are already inhibited by the formation of the PD-L1-PD-1 complex. This complex forms due to T cell interaction with tumor cells, and it masks the tumor from identification by T cells. Incorporating these two forms of immunosuppression prepares our model for future extension to include immunotherapies and optimize treatments.
Our results reinforce pre-clinical studies by Flores-Toro et al. (2020) which suggest that MDSCs and the PD-L1-PD-1 immune checkpoint should be targeted together to increase the immune response and thus decrease the tumor burden. Therefore, promising therapeutics include immune checkpoint inhibitors, such as anti-PD-1 and anti-PD-L1, and CCR2 antagonists. Time-dependent sensitivity analysis with eFAST indicates that treatment with these therapeutics should be used as soon as possible; however, an immune checkpoint inhibitor should offer some benefit at any time. In contrast, MDSC recruitment should be blocked by a CCR2 antagonist in the early stages of a tumor, but immunotherapies which kill MDSCs at the tumor site should be considered for larger tumors. While CCR2 antagonists prevent MDSC recruitment (thus targeting the \(s_M\) parameter), strategies which directly target the parameter r by decreasing the ability of MDSCs to inhibit T cells should be explored. Time-dependent eFAST results suggest that therapies targeting r would be especially useful during the later stages of tumor development.
Stability analysis of the system (1) established conditions for the existence and stability of the tumor-free equilibrium, \((0,T^*_0, 0)\). This unique equilibrium is corroborated by immunohistochemistry data showing the absence of MDSCs within the tumor-free mouse brain (Fig. 3). When \(\lambda _C < \eta T_0^*\), i.e., the initial immune response is greater than the tumor growth rate, the tumor-free equilibrium is locally asymptotically stable, which is equivalent to results in Nikolopoulou et al. (2018). Further, it is globally stable when \(\lambda _C (1+\rho \widehat{T}(\widehat{T} +\epsilon _C \widehat{C}))(r\widehat{M} +d_T) < \eta a_T \). This suggests that, regardless of the initial number of activated T cells, if T cells are activated and kill tumor cells faster than tumor cells proliferate and faster than T cells are suppressed and die naturally, then the system will tend to the tumor-free equilibrium.
Bifurcation analysis further illuminates this result by showing that shifting only 4 parameters (\(\lambda _C\), \(\eta \), \(\rho \), r) can result in the system tending to the tumor-free equilibrium. As surgical resection of a tumor can decrease the tumor growth rate by disrupting the influx of nutrients through damage to the vasculature, these results suggest that a treatment plan including surgery and immunotherapies which target the PD-L1-PD-1 complex and the immunosuppressive capabilities of MDSCs would minimize the tumor size.
We obtained estimates of the tumor volume and cell counts from murine GBM data collected from immunofluorescence imaging. T cell data was more sparse than the glioma and MDSC data, so future model predictions could be improved by incorporating more time-dependent T cell data. Numerical simulations in Fig. 5a show that our parameter distributions capture nearly all data points (33 of 36 data points) within one standard deviation from the mean and most (27 of 36) within half of a standard derivation. In general, more data would improve the predictive capabilities of the model, especially if the cell populations within a single mouse could be tracked over time.
Given the limited data set, we used the Approximate Bayesian Computation (ABC) rejection method in unison with data to obtain information on the probability distribution of realistic parameter values. We identified approximate values and quantified the associated uncertainty for previously unknown parameters, in particular, the recruitment rate of MDSCs (\(s_M\)) and the inhibition rate of T cells by MDSCs (r). By calculating an estimate of \(s_M\) using chemokine expression levels on glioma cells and MDSC migration in response to certain concentrations of the chemokine CCL2 measured from experimental data, we found that the parameter values estimated by the ABC method were of the same order as calculated from data. Further, we found that the values estimated for r were noticeably greater than the estimated values for \(\eta \) (Table 2), suggesting that MDSCs suppress T cells more effectively than T cells kill tumor cells. This would explain the difficulty T cells experience in overcoming a glioma. Further biological testing would need to occur to validate this hypothesis.
The ABC parameter distributions (Fig. 4) informed global sensitivity analysis by identifying a search curve with which to sample parameters for eFAST. Sensitivity analysis (Fig. 6) indicated that there are several drivers of the system, including the tumor growth rate (\(\lambda _C\)), tumor carrying capacity (\(C_\textrm{max}\)), and T cell kill rate (\(\eta \)), along with noticeable interaction occurring between parameters to affect the variance of the system. By comparing the eFAST results with the ABC results, we hypothesize that interaction could be occurring between these three drivers and parameters representing immune suppression–namely the tumor upregulation of PD-L1 (\(\epsilon _C\)) along with the inhibition by PD-L1-PD-1 (\(\rho \)) and by MDSCs (r), since there appears to be an inverse relationship between these parameters in Fig. 4a. Further, T cell sensitivity in Fig. 6 shows that targeting these parameters would increase the immune response regardless of any interactions. This suggests that the PD-L1-PD-1 complex and inhibition by MDSCs should be therapeutically targeted together to relieve tumor burden and increase the immune response.
Overall, there is a mutually beneficial relationship between the ABC method and eFAST. The distributions produced by ABC in accordance with data inform the choice of search curve when sampling parameters in eFAST. In return, the interaction effects displayed by eFAST direct us to look for relationships between parameters using the ABC results. Thus, the two together increase the reliability of results as well as the inferences that can be made.
Future directions include extending this model to incorporate immunotherapies targeting the PD-L1-PD-1 complex and MDSCs and then optimizing treatment regimens to minimize tumor burden. Anti-PD-1, which targets PD-1 on T cells, has not been successful as a monotherapy in part because it requires more activated T cells to be in the tumor microenvironment (Kleponis et al. 2015). Simulations in Nikolopoulou et al. (2018) suggest that anti-PD-1 alone is not able to increase the number of T cells to a tumor eradication threshold. Through inhibition, MDSCs decrease the number of activated T cells. Treatment with a CCR2 antagonist decreases the recruitment of MDSCs, which indirectly increases the number of activated T cells, thus improving the efficacy of anti-PD-1. This improved efficacy is supported by data showing increased survival rates in glioma-bearing mice treated with anti-PD-1 and the CCR2 antagonist, CCX872 (Flores-Toro et al. 2020, Fig. 4). While our results agree with this outcome, we also propose an additional target. Collectively, our global stability and bifurcation analyses, sensitivity analysis, and ABC method results suggest that it would be more beneficial to directly target the mechanisms by which MDSCs inhibit T cells rather than MDSC recruitment. Therefore, therapeutics which prevent MDSC inhibition of T cells should be tested in combination with immune checkpoint inhibitors like anti-PD-1.
Code availability
The source code used to generate the results for this article is available through GitHub at https://github.com/stepien-lab/glioma-Tcell-MDSC [v1.0.0]. The code is platform independent and written in MATLAB.
Change history
References
Allahverdy A, Moghaddam AK, Rahbar S et al (2019) An agent-based model for investigating the effect of myeloid-derived suppressor cells and its depletion on tumor immune surveillance. J Med Signals Sens 9(1):15–23. https://doi.org/10.4103/jmss.JMSS_33_18
Allmaras M, Bangerth W, Linhart JM et al (2013) Estimating parameters in physical models through Bayesian inversion: a complete example. Soc Ind Appl Math Rev 55(1):149–167. https://doi.org/10.1137/100788604
Ansell SM, Lesokhin AM, Borrello I et al (2015) PD-1 blockade with Nivolumab in relapsed or refractory Hodgkin’s lymphoma. N Engl J Med 372(4):311–319. https://doi.org/10.1056/NEJMoa1411087
Banerjee S, Khajanchi S, Chaudhuri S (2015) A mathematical model to elucidate brain tumor abrogation by immunotherapy with T11 target structure. PLoS ONE 10(5):e0123611. https://doi.org/10.1371/journal.pone.0123611
Bausart M, Préat V, Malfanti A (2022) Immunotherapy for glioblastoma: the promise of combination strategies. J Exp Clinical Cancer Res 41(1):1–22. https://doi.org/10.1186/s13046-022-02251-2
Bitsouni V, Tsilidis V (2022) Mathematical modeling of tumor-immune system interactions: the effect of rituximab on breast cancer immune response. J Theor Biol 539(111):001. https://doi.org/10.1016/j.jtbi.2021.111001
Brahmer J, Reckamp KL, Baas P et al (2015) Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer. N Engl J Med 373(2):123–135. https://doi.org/10.1056/NEJMoa1504627
Brown NF, Carter TJ, Ottaviani D et al (2018) Harnessing the immune system in glioblastoma. Br J Cancer 119(10):1171–1181. https://doi.org/10.1038/s41416-018-0258-8
Browning AP, McCue SW, Simpson MJ (2017) A Bayesian computational approach to explore the optimal duration of a cell proliferation assay. Bull Math Biol 79(8):1888–1906. https://doi.org/10.1007/s11538-017-0311-4
Bryukhovetskiy I (2022) Cell-based immunotherapy of glioblastoma multiforme. Oncol Lett 23(4):1–14. https://doi.org/10.3892/ol.2022.13253
Butner JD, Wang Z, Elganainy D et al (2021) A mathematical model for the quantification of a patient’s sensitivity to checkpoint inhibitors and long-term Tumour burden. Nat Biomed Eng 5(4):297–308. https://doi.org/10.1038/s41551-020-00662-0
Cao Y, Feng Y, Zhang Y et al (2016) 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 16(1):1–11. https://doi.org/10.1186/s12885-016-2376-0
Chang AL, Miska J, Wainwright DA et al (2016) CCL2 produced by the glioma microenvironment is essential for the recruitment of regulatory T cells and myeloid-derived suppressor cells. Can Res 76(19):5671–5682. https://doi.org/10.1158/0008-5472.CAN-16-0144
Chung H, Ros W, Delord J et al (2019) Efficacy and safety of pembrolizumab in previously treated advanced cervical cancer: results from the phase II KEYNOTE-158 study. J Clin Oncol 37(17):1470–1478. https://doi.org/10.1200/JCO.18.01265
Condamine T, Gabrilovich DI (2011) Molecular mechanisms regulating myeloid-derived suppressor cell differentiation and function. Trends Immunol 32(1):19–25. https://doi.org/10.1016/j.it.2010.10.002
da Costa JMJ, Orlande HRB, da Silva WB (2018) Model selection and parameter estimation in tumor growth models using approximate Bayesian computation-ABC. Comput Appl Math 37(3):2795–2815. https://doi.org/10.1007/S40314-017-0479-0
Cukier R, Fortuin C, Shuler KE et al (1973) Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients I Theory. J Chem Phys 59(8):3873–3878. https://doi.org/10.1063/1.1680571
Doblas S, He T, Saunders D et al (2010) Glioma morphology and tumor-induced vascular alterations revealed in seven rodent glioma models by in vivo magnetic resonance imaging and angiography. J Magn Reson Imaging 32(2):267–275. https://doi.org/10.1002/jmri.22263
Duraiswamy J, Kaluza KM, Freeman GJ et al (2013) Dual blockade of PD-1 and CTLA-4 combined with tumor vaccine effectively restores T-cell rejection function in tumors. Can Res 73(12):3591–3603. https://doi.org/10.1158/0008-5472.CAN-12-4100
Eftimie R, Bramson JL, Earn DJ (2011) Interactions between the immune system and cancer: a brief review of non-spatial mathematical models. Bull Math Biol 73(1):2–32. https://doi.org/10.1007/s11538-010-9526-3
Fernandes C, Costa A, Osório L, et al (2017) Glioblastoma [Internet], Codon Publications, Brisbane, Australia, chap 11: Current standards of care in glioblastoma therapy, pp 197–241. https://doi.org/10.15586/codon.glioblastoma.2017.ch11
Flores-Toro JA, Luo D, Gopinath A et al (2020) CCR2 inhibition reduces tumor myeloid cells and unmasks a checkpoint inhibitor effect to slow progression of resistant murine gliomas. Proc Natl Acad Sci 117(2):1129–1138. https://doi.org/10.1073/pnas.1910856117
Gabrilovich DI, Nagaraj S (2009) Myeloid-derived suppressor cells as regulators of the immune system. Nat Rev Immunol 9(3):162–174. https://doi.org/10.1038/nri2506
Grover WH, Bryan AK, Diez-Silva M et al (2011) Measuring single-cell density. Proc Natl Acad Sci 108(27):10992–10996. https://doi.org/10.1073/pnas.1104651108
Jafarnejad M, Gong C, Gabrielson E et al (2019) A computational model of neoadjuvant PD-1 inhibition in non-small cell lung cancer. J Amer Assoc Pharm Sci 21(5):1–14. https://doi.org/10.1208/s12248-019-0350-x
Khajanchi S (2021) The impact of immunotherapy on a glioma immune interaction model. Chaos, Solitons Fractals 152(111):346. https://doi.org/10.1016/j.chaos.2021.111346
Khajanchi S, Banerjee S (2017) Quantifying the role of immunotherapeutic drug T11 target structure in progression of malignant gliomas: mathematical modeling and dynamical perspective. Math Biosci 289:69–77. https://doi.org/10.1016/j.mbs.2017.04.006
Khyat T, Jang SRJ (2022) On a discrete model of Tumour-immune system interactions with blockade of immune checkpoints. J Differ Equ Appl 28(1):73–108. https://doi.org/10.1080/10236198.2021.2023136
Kirschner D (2008) Uncertainty and sensitivity functions and implementation. http://malthus.micro.med.umich.edu/lab/usadata/
Kirschner D, Panetta JC (1998) Modeling immunotherapy of the tumor-immune interaction. J Math Biol 37(3):235–252. https://doi.org/10.1007/s002850050127
Kleponis J, Skelton R, Zheng L (2015) Fueling the engine and releasing the break: combinational therapy of cancer vaccines and immune checkpoint inhibitors. Cancer Biol Med 12(3):201–208. https://doi.org/10.7497/j.issn.2095-3941.2015.0046
Korbecki J, Kojder K, Simińska D et al (2020) CC chemokines in a tumor: a review of pro-cancer and anti-cancer properties of the ligands of receptors CCR1, CCR2, CCR3, and CCR4. Int J Mol Sci 21(21):8412. https://doi.org/10.3390/ijms21218412
Kreger J, Roussos Torres ET, MacLean AL (2023) Myeloid-derived suppressor-cell dynamics control outcomes in the metastatic niche. Cancer Immunol Res 11(5):614–628. https://doi.org/10.1158/2326-6066.CIR-22-0617
Kuznetsov VA, Makalkin IA, Taylor MA et al (1994) Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bull Math Biol 56(2):295–321. https://doi.org/10.1016/S0092-8240(05)80260-5
Lai X, Friedman A (2017) Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitors: A mathematical model. PLoS ONE 12(5):e0178,479. https://doi.org/10.1371/journal.pone.0178479
Lai X, Stiff A, Duggan M et al (2018) Modeling combination therapy for breast cancer with BET and immune checkpoint inhibitors. Proc Natl Acad Sci 115(21):5534–5539. https://doi.org/10.1073/pnas.1721559115
Liao KL, Bai XF, Friedman A (2014) Mathematical modeling of interleukin-35 promoting tumor growth and angiogenesis. PLoS ONE 9(10):e110,126. https://doi.org/10.1371/journal.pone.0110126
Liepe J, Kirk P, Filippi S et al (2014) A framework for parameter estimation and model selection from experimental data in systems biology using approximate Bayesian computation. Nat Protoc 9(2):439–456. https://doi.org/10.1038/nprot.2014.025
Lim M, Xia Y, Bettegowda C et al (2018) Current state of immunotherapy for glioblastoma. Nat Rev Clin Oncol 15(7):422–442. https://doi.org/10.1038/s41571-018-0003-5
Liu C, Yu S, Kappes J et al (2007) Expansion of spleen myeloid suppressor cells represses NK cell cytotoxicity in tumor-bearing host. Blood 109(10):4336–4342. https://doi.org/10.1182/blood-2006-09-046201
Mahlbacher GE, Reihmer KC, Frieboes HB (2019) Mathematical modeling of tumor-immune cell interactions. J Theor Biol 469:47–60. https://doi.org/10.1016/j.jtbi.2019.03.002
Markowitz J, Wang J, Vangundy Z et al (2017) Nitric oxide mediated inhibition of antigen presentation from DCs to CD4+ T cells in cancer and measurement of STAT1 nitration. Sci Rep 7(1):1–13. https://doi.org/10.1038/s41598-017-14970-0
Mirzaei N, Su S, Sofia D et al (2021) A mathematical model of breast tumor progression based on immune infiltration. J Personalized Med 11(10):1031. https://doi.org/10.3390/jpm11101031
Monu NR, Frey AB (2012) Myeloid-derived suppressor cells and anti-tumor T cells: a complex relationship. Immunol Invest 41(6–7):595–613. https://doi.org/10.3109/08820139.2012.673191
Nikolopoulou E, Johnson L, Harris D et al (2018) Tumour-immune dynamics with an immune checkpoint inhibitor. Lett Biomath 5(2):S137–S159. https://doi.org/10.30707/LiB5.2Nikolopoulou
Nikolopoulou E, Eikenberry SE, Gevertz JL et al (2021) Mathematical modeling of an immune checkpoint inhibitor and its synergy with an immunostimulant. Discr Contin Dyn Syst B 26(4):2133. https://doi.org/10.3934/dcdsb.2020138
Ostrum Q, Gittleman H, Fulop J et al (2015) CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2008–2012. Neuro Oncol 17(Suppl 4):iv1–iv62. https://doi.org/10.1093/neuonc/nov189
Özköse F, Yılmaz S, Yavuz M et al (2022) A fractional modeling of tumor-immune system interaction related to lung cancer with real data. Euro Phys J Plus 137(1):1–28. https://doi.org/10.1140/epjp/s13360-021-02254-6
Perlstein D, Shlagman O, Kogan Y et al (2019) Personal response to immune checkpoint inhibitors of patients with advanced melanoma explained by a computational model of cellular immunity, tumor growth, and drug. PLoS ONE 14(12):e0226,869. https://doi.org/10.1371/journal.pone.0226869
Pillay J, Den Braber I, Vrisekoop N et al (2010) In vivo labeling with 2H2O reveals a human neutrophil lifespan of 5.4 days. Blood 116(4):625–627. https://doi.org/10.1182/blood-2010-01-259028
Postow MA, Chesney J, Pavlick AC et al (2015) Nivolumab and Ipilimumab versus Ipilimumab in untreated melanoma. N Engl J Med 372(21):2006–2017. https://doi.org/10.1056/NEJMoa1414428
Preusser M, Lim M, Hafler DA et al (2015) Prospects of immune checkpoint modulators in the treatment of glioblastoma. Nat Rev Neurol 11(9):504–514. https://doi.org/10.1038/nrneurol.2015.139
Radunskaya A, Kim R, Woods T II (2018) Mathematical modeling of tumor immune interactions: a closer look at the role of a PD-L1 inhibitor in cancer immunotherapy. SPORA: A J Biomath 4(1):25–41. https://doi.org/10.30707/SPORA4.1Radunskaya
Reardon D, Omuro A, Brandes A et al (2017) OS10. 3 randomized phase 3 study evaluating the efficacy and safety of Nivolumab vs bevacizumab in patients with recurrent glioblastoma: Checkmate 143. Neuro Oncol 19(Suppl 3):iii21. https://doi.org/10.1093/neuonc/nox036.071
Ribeiro RM, Mohri H, Ho DD et al (2002) In vivo dynamics of T cell activation, proliferation, and death in HIV-1 infection: why are CD4+ but not CD8+ T cells depleted? Proc Natl Acad Sci 99(24):15572–15577. https://doi.org/10.1073/pnas.242358099
Rutter EM, Stepien TL, Anderies BJ et al (2017) Mathematical analysis of glioma growth in a murine model. Sci Rep 7(2508):1–16. https://doi.org/10.1038/s41598-017-02462-0
Saio M, Radoja S, Marino M et al (2001) Tumor-infiltrating macrophages induce apoptosis in activated CD8+ T cells by a mechanism requiring cell contact and mediated by both the cell-associated form of TNF and nitric oxide. J Immunol 167(10):5583–5593. https://doi.org/10.4049/jimmunol.167.10.5583
Saltelli A, Tarantola S, Chan KS (1999) A quantitative model-independent method for global sensitivity analysis of model output. Technometrics 41(1):39–56. https://doi.org/10.1080/00401706.1999.10485594
Saltelli A, Ratto M, Andres T et al (2008) Global sensitivity analysis: the primer. John Wiley & Sons Ltd, Chichester, England. https://doi.org/10.1002/9780470725184
Shariatpanahi SP, Shariatpanahi SP, Madjidzadeh K et al (2018) Mathematical modeling of tumor-induced immunosuppression by myeloid-derived suppressor cells: Implications for therapeutic targeting strategies. J Theor Biol 442:1–10. https://doi.org/10.1016/j.jtbi.2018.01.006
Shi S, Huang J, Kuang Y (2021) Global dynamics in a tumor-immune model with an immune checkpoint inhibitor. Discr Contin Dyn Syst B 26(2):1149–1170. https://doi.org/10.3934/dcdsb.2020157
Srivastava MK, Zhu L, Harris-White M et al (2012) Myeloid suppressor cell depletion augments antitumor activity in lung cancer. PLoS ONE 7(7):e40677. https://doi.org/10.1371/journal.pone.0040677
Stensjøen AL, Solheim O, Kvistad KA et al (2015) Growth dynamics of untreated glioblastomas in vivo. Neuro Oncol 17(10):1402–1411. https://doi.org/10.1093/neuonc/nov029
Stepien TL, Lynch HE, Yancey SX et al (2019) Using a continuum model to decipher the mechanics of embryonic tissue spreading from time-lapse image sequences: An approximate Bayesian computation approach. PLoS ONE 14(6):e0218021. https://doi.org/10.1371/journal.pone.0218021
Storey KM, Lawler SE, Jackson TL (2020) Modeling oncolytic viral therapy, immune checkpoint inhibition, and the complex dynamics of innate and adaptive immunity in glioblastoma treatment. Front Physiol 11(151):00151. https://doi.org/10.3389/fphys.2020.00151
Stupp R, Mason WP, Van Den Bent MJ et al (2005) Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med 352(10):987–996. https://doi.org/10.1056/NEJMoa043330
Sunnåker M, Busetto AG, Numminen E et al (2013) Approximate Bayesian computation. PLoS Comput Biol 9(1):e1002803. https://doi.org/10.1371/journal.pcbi.1002803
Takacs GP, Flores-Toro JA, Harrison JK (2021) Modulation of the chemokine/chemokine receptor axis as a novel approach for glioma therapy. Pharmacol Therapeut 222(107):790. https://doi.org/10.1016/j.pharmthera.2020.107790
Takacs GP, Kreiger CJ, Luo D et al (2022) Glioma-derived CCL2 and CCL7 mediate migration of immune suppressive CCR2+/CX3CR1+ M-MDSCs into the tumor microenvironment in a redundant manner. Front Immunol 13:7959. https://doi.org/10.3389/fimmu.2022.993444
Vetsika EK, Koukos A, Kotsakis A (2019) Myeloid-derived suppressor cells: major figures that shape the immunosuppressive and angiogenic network in cancer. Cells 8(12):1647. https://doi.org/10.3390/cells8121647
Xiao Y, Thomas L, Chaplain MA (2021) Calibrating models of cancer invasion: parameter estimation using approximate Bayesian computation and gradient matching. Royal Soc Open Sci 8(202):237. https://doi.org/10.1098/rsos.202237
Yu JL, Jang SRJ (2019) A mathematical model of tumor-immune interactions with an immune checkpoint inhibitor. Appl Math Comput 362(C):1–11. https://doi.org/10.1016/j.amc.2019.06.037
Yu MW, Quail DF (2021) Immunotherapy for glioblastoma: current progress and challenges. Front Immunol 12(676):301. https://doi.org/10.3389/fimmu.2021.676301
Acknowledgements
The authors would like to thank and acknowledge the reviewers for the time, effort, expertise, and comments given on the original manuscript.
Funding
H.G.A. and G.P.T. acknowledge support from the NIH/NCATS Clinical and Translational Science Awards #UL1TR001427 and #TL1TR001428. H.G.A. acknowledges support from NSF grant DMS-2151566. D.C.H. acknowledges support from NIH R01GM131405-02. Y.K. acknowledges support from NSF DEB-1930728 and NIH R01GM131405-02. J.K.H. acknowledges support from the National Institute of Neurological Disorders and Stroke NS108781. T.L.S. acknowledges support from a Simons Collaboration Grant for Mathematicians (#710482) and NSF grant DMS-2151566.
Author information
Authors and Affiliations
Contributions
Conceptualization and Methodology: all authors; Formal analysis, Software, Validation, and Visualization: Hannah G. Anderson; Investigation and Data Curation: Gregory P. Takacs; Resources: Jeffrey K. Harrison; Supervision: Tracy L. Stepien, Jeffrey K. Harrison, Yang Kuang; Writing—original draft preparation: Hannah G. Anderson; Writing—review and editing: all authors.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Ethical approval
All procedures involving animal housing and surgical protocols were followed according to the guidelines of the University of Florida Institutional Animal Care and Use Committee.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Below is the link to the electronic supplementary material.
Online Resource 1:
ESM1.pdf contains details regarding data collection and usage, the literature search and estimation of parameters, as well as additional figures and information about the ABC rejection method, bifurcation diagrams, and eFAST. (pdf 422KB)
Online Resource 2:
ESM2.xlsx contains time-dependent experimental data of tumor volumes and cell counts of the three studied cell populations, which were collected according to the methods described in Online Resource 1. (xl 18KB)
Appendix A Additional Theorems
Appendix A Additional Theorems
In the following appendix, we establish conditions for the existence and stability of a unique tumorous equilibrium.
Theorem 6
(Uniqueness of tumorous equilibrium) The system (1) has a tumorous equilibrium if \(\lambda _C > \eta T^*\) and \(C_\textrm{max} \epsilon _C \eta /\lambda _C <1\). Further, it is unique when \(\epsilon _C \rho > \eta /\lambda _C\), \(C_\textrm{max} \ge 2\), \(\beta _1/\beta _2 < C_\textrm{max} \epsilon _C \lambda _C/\eta \), and \(\rho \beta _1 \left( \lambda _C/\eta \right) ^2 < C_\textrm{max} \big (s_M r - d_M s_T\big )\), where \(\beta _i=d_Td_M + \alpha r +s_M r\big ((i+1)C_\textrm{max} +q\big )\), for \(i=1,2\).
Proof
Setting the right hand side of the tumor cell Eq. (1a) equal to zero and evaluating at the tumorous equilibrium implies that
which is positive if
For the MDSCs Eq. (1c), we find that
which is positive if \(C^*\) is positive.
The T cell Eq. (1b) yields
which, when we substitute (A3) and (A1) and rearrange, becomes a degree-5 polynomial,
where
We observe that \(a_0\) is always positive since all the parameters are positive. Since \(a_0\) is positive, when we require \(a_5\) to be negative, the number of sign changes is \(1 \ \textrm{mod} \ 2\) regardless of the signs of the coefficients \(a_1,\ldots ,a_4\). Therefore, by Descartes’ rule of signs, there is at least one positive zero for (A5). Thus, a tumorous equilibrium is guaranteed to exist if we enforce condition (A2) and
We derive conditions for uniqueness in the case where \(a_i<0\) for \(i=1,\ldots ,5\). \(a_1\) is negative if we require
In other words, immunosuppression due to MDSCs must be greater than the recruitment of T cells and death of MDSCs.
For \(a_2\) to be negative, under the conditions that
it follows that
When we require
then
We shall consider \(a_4\) next. When we expand \(a_4\) by distributing the term containing \(d_Td_M\), it becomes clearer to see that \(a_4\) is negative under the condition that
For \(a_3\), conditions (A8) and (A7) cause the first and second terms of \(a_3\) to be negative, respectively, while (A13) implies that the last term is positive. Applying (A7) to the last three terms and simplifying,
Thus, \(a_3\) is negative under the condition that
which is a stronger version of (A8).
Under this condition along with (A7), (A9), (A11), and (A13), we have \(a_i < 0\) for \(i = 1,...,5\) and \(a_0 >0\). Thus, according to Descartes’ rule of signs, there is one unique positive zero for \(T^*\), and correspondingly from (A1) and (A3), one unique \(C^*\) and \(M^*\). Therefore, we have a unique tumorous equilibrium \((C^*, T^*, M^*)\). \(\square \)
There are some biological implications from these conditions. First of all, condition (A9) can be rewritten as \(\eta < \lambda _C \epsilon _C \rho \). This suggests that the tumor growth rate (\(\lambda _C\)) multiplied by immunosuppression via the PD-L1-PD-1 complex (\(\epsilon _C \rho \)) must be larger than the ability of T cells to kill tumor cells (\(\eta \)). Further, condition (A15) indicates that \(d_M s_T < s_M r\), in other words, immunosuppression due to MDSCs (\(s_M r\)) must be greater than the recruitment of T cells (\(s_T\)) and death of MDSCs (\(d_M\)). Both implications suggest that the suppression of the immune system has to be substantial enough for uniqueness of the tumorous equilibrium to be guaranteed.
Next we determine conditions for the stability of the tumorous equilibrium.
Theorem 7
(Stability of tumorous equilibrium) If a tumorous equilibrium, \((C,T,M)=(C^*, T^*, M^*)\), exists as determined by the conditions in Theorem 6, it is locally asymptotically stable when \(\lambda _C > 2 \eta T^*\) and \(\rho T^{*2} \ge 1\).
Proof
The Jacobian evaluated at the tumorous equilibrium is
where
and we have substituted (A1) into the (1, 1)-entry to simplify the expression.
When we require
k is guaranteed to be negative.
The characteristic polynomial of \(J(C^*,T^*,M^*)\) is
A monic polynomial is defined to be Hurvitz if \(\text {Re}(x_i)<0\) for all roots \(x_i \in \mathbb {C}\). By the Routh-Hurvitz Criterion, a polynomial \(p(x) = x^3 +b_1 x^2 + b_2 x + b_3\) is Hurvitz if and only if \(b_i >0\) for \(i=1,2,3\), and \((b_1 b_2 -b_3) >0\).
Under conditions (A2) and (A18) it follows that
If we strengthen condition (A2) to be
this allows us to conclude, after substituting g and k and simplifying, that
Therefore,
Next, since \(rM^*<-k\), and substituting in \(M^*\) (A3), it follows by condition (A21) that
Finally, we shall address the positivity of \((b_1b_2 -b_3)\),
Now, \(\lambda _C -\eta T^* -k>0\) since \(k<0\) and by condition (A21). Therefore, by condition (A22),
Since \(b_i >0\) for \(i=1,2,3\) and \((b_1 b_2 -b_3) >0\), the characteristic polynomial (A19) is Hurvitz. Therefore, the tumorous equilibrium is locally asymptotically stable under the conditions (A18) and (A21). \(\square \)
Condition (A18) requires that the T cell inhibition by the PD-L1-PD-1 complex is at or above a certain level. Condition (A21) effectively assumes that the tumor cell growth rate is larger than twice the kill rate by the T cell population. We note that Nikolopoulou et al. (2018) also required a similar condition (\(\lambda _C > \eta T^*\)) in order to achieve a stable tumorous equilibrium for their system.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Anderson, H.G., Takacs, G.P., Harris, D.C. et al. Global stability and parameter analysis reinforce therapeutic targets of PD-L1-PD-1 and MDSCs for glioblastoma. J. Math. Biol. 88, 10 (2024). https://doi.org/10.1007/s00285-023-02027-y
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00285-023-02027-y
Keywords
- Immune checkpoint
- Immunosuppression
- Mathematical oncology
- Approximate Bayesian computation
- Sensitivity analysis