Abstract
The emergence of drug resistance is often an inevitable obstacle that limits the long-term effectiveness of clinical cancer chemotherapeutics. Although various forms of cancer cell-intrinsic mechanisms of drug resistance have been experimentally revealed, the role and the underlying mechanism of tumor microenvironment in driving the development of acquired drug resistance remain elusive, which significantly impedes effective clinical cancer treatment. Recent experimental studies have revealed a macrophage-mediated drug resistance mechanism in which the tumor microenvironment undergoes adaptation in response to macrophage-targeted colony-stimulating factor-1 receptor (CSF1R) inhibition therapy in gliomas. In this study, we developed a spatio-temporal model to quantitatively describe the interplay between glioma cells and CSF1R inhibitor–targeted macrophages through CSF1 and IGF1 pathways. Our model was used to investigate the evolutionary kinetics of the tumor regrowth and the associated dynamic adaptation of the tumor microenvironment in response to the CSF1R inhibitor treatment. The simulation result obtained using this model was in agreement with the experimental data. The sensitivity analysis revealed the key parameters involved in the model, and their potential impacts on the model behavior were examined. Moreover, we demonstrated that the drug resistance is dose-dependent. In addition, we quantitatively evaluated the effects of combined CSFR inhibition and IGF1 receptor (IGF1R) inhibition with the goal of designing more effective therapies for gliomas. Our study provides quantitative and mechanistic insights into the microenvironmental adaptation mechanisms that operate during macrophage-targeted immunotherapy and has implications for drug dose optimization and the design of more effective combination therapies. Mol Cancer Ther; 17(4); 814–24. ©2018 AACR.
Introduction
As a microenvironment-targeted immunotherapy, colony-stimulating factor-1 receptor (CSF1R) inhibition is expected to be a promising therapeutic strategy for high-grade gliomas (1, 2). Although multiple clinical trials testing the efficacy of CSF1R inhibition have been conducted in glioma patients, it is unclear whether and how resistance emerges during this therapy. Recently, a preclinical study using a genetic mouse model (3) showed that a high percentage (>50%) of glioblastoma multiformes (GBM) recur in the mice following treatment via long-term CSF1R inhibition, although overall survival is significantly prolonged. Therefore, it is important and timely to investigate the mechanisms underlying the acquisition of resistance to CSF1R inhibition and quantitatively determine the key biological parameters involved in this process.
Many types of tumor cell-intrinsic mechanisms, including genetic mutations and epigenetic modifications (4), signaling network rewiring (5–7), cellular heterogeneity involving cancer stem cells (8, 9), and drug metabolism and pharmacodynamics (10), have been studied widely. Distinct from the above mechanisms, the roles of tumor microenvironment in driving and facilitating tumor growth and drug resistance have been revealed and acknowledged with increasing attention in recent years (3, 11–13). Actually, the tumor microenvironment involves multiple cell types, such as fibroblasts, endothelial cells and immune cells, that can interact with tumor cells (14). For example, immune cells, including macrophages, dendritic cells, neutrophils, Treg cells, T effector cells, natural killer cells, etc. (15), can interact with tumor cells either via direct cell–cell contacts or via soluble factors such as secreted cytokines or growth factors (16). Such cross-talk between immune and tumor cells plays substantial roles before and during tumor growth as well as the in the development of drug resistance (17).
To investigate the dynamic interactions between cancer cells and the associated microenvironment in the context of therapeutic intervention, the traditional experimental approach requires multiple experimental conditions and time points and, thus, is time consuming and expensive. Therefore, a computational systems biology approach is necessary to address this issue. Mathematical modeling is a powerful tool for the quantitative investigation of dynamic biological systems that can provide insights into underlying mechanisms and experimentally testable predictions. From a quantitative biology perspective, we can ask questions regarding the effects of dose on drug resistance and the influence of critical parameters on the evolution of drug resistance. Driven by the hypothesis that an appropriately simplified mathematical model is able to capture the key kinetics of microenvironment-mediated drug resistance, we seek to mechanistically model the spatio-temporal evolution of interactions within the tumor-microenvironment system during drug treatment.
In the past decade, several theoretical models have been developed to investigate the role of the microenvironment in tumor growth [e.g., refs. (18–21)]. Moreover, some efforts using modeling approaches to investigate the role of the microenvironment in drug resistance have also been conducted recently. For instance, Lindsay and colleagues (22) developed a stochastic model to simulate the effect of oxygen on the growth kinetics of cancer cells and to investigate the potential benefits of combining hypoxia-activated prodrugs with standard targeted therapy to prevent drug resistance in non–small cell lung cancer. In addition, our previous study (23) designed a modeling framework using stochastic differential equations to simulate the therapy-induced microenvironmental adaptation that promotes the growth and metastasis of resistant tumor cells, and therefore contributes to the development of the acquired drug resistance. However, no modeling work has been dedicated to investigating the kinetics of drug resistance driven by the immune microenvironment, particularly involving macrophages, in cancer therapy.
Recently, the experimental studies have revealed a macrophage-mediated drug resistance mechanism in glioma immunotherapy (3). Specifically, in the absence of CSF1R inhibition, tumor-associated macrophages exhibit a protumorigenic M2 phenotype and contribute to glioma progression. However, when targeted by CSF1R inhibitors, macrophages switch to an antitumorigenic M1 phenotype with enhanced phagocytosis, thereby promoting tumor cell death. Following prolonged treatment with CSF1R inhibitors, accumulated IL4 from other cell types (e.g., T cells) stimulates macrophages to secrete insulin-like growth factor 1 (IGF1) into the microenvironment, which in turn sustains the survival and growth of glioma cells (3).
On the other hand, certain mathematical models have also recently been developed to describe the interactions between cancer cells and macrophages during cancer progression. For example, Szomolay and colleagues (24) proposed a system of partial differential equations (PDE) to simulate the macrophages-mediated interactions among tumor cells and endothelial cells through several cytokines, including M-CSF1, GM-CSF1, VEGF, and MCP-1, as well as oxygen. Their model simulations are consistent with the in vivo data obtained under different conditions of GM-CSF (granulocyte/macrophage colony-stimulating factor) exposure, allowing the model to predict treatment outcomes. Knútsdóttir and colleagues (25) described paracrine and autocrine signaling (i.e., CSF1/CSF1R and EGF/EGFR pathways) between tumor cells and macrophages using both PDEs model and cell-based discrete 3D simulation. The authors identified parameters responsible for the observed ratio of tumor cell to macrophage and demonstrated how CSF1/CSF1R autocrine signaling in tumor cells affect this ratio. These works provided the modeling basis and implicated potential values of theoretically investigating the interactions between tumor cells and macrophages.
In this study, based on the above experimentally validated mechanism (3), we developed a spatio-temporal model to quantitatively describe the interplay between glioma cells and CSF1R inhibition-induced transition of macrophage M2 and M1 phenotypes through dynamic CSF1/CSF1R and IGF1/IGF1R pathways. Our mode was used to investigate the response of glioma cells to macrophage-targeted CSF1R inhibition therapy and the associated dynamic adaptation of the tumor microenvironment. The simulation result obtained using this model is in agreement with the experimentally observed regrowth phenomenon of resistant gliomas and quantitative survival times (3). Moreover, we demonstrated that the drug resistance is dose-dependent. In addition, we quantitatively evaluated the effects of combined CSFR inhibition and IGF1 receptor (IGF1R) inhibition with the goal of designing more effective therapies for gliomas.
Materials and Methods
PDE model
We model the spatio-temporal evolution of cancer cells and their microenvironment during CSF1R inhibition based on the experimentally validated mechanism (ref. 3; Fig. 1). The major assumptions of our model include:
a) The numbers of cancer cells and macrophages as well as levels of secreted soluble factors including CSF1, IGF1, and IL4 are large enough such that they can be treated as continuous variables. So, the numbers of cells and levels of cytokines are described by their densities and concentrations, respectively.
b) Glioma cells, macrophages, CSF1, IGF1, and drugs are subject to diffusion; glioma cells and macrophages are also subject to chemotaxis in response to IGF1 and CSF1, respectively.
c) The proliferation rate and the death rate of tumor cells are regulated by IGF1 and CSF1R inhibitor. Specifically, by regulating the phenotype transition between M1 and M2 macrophages, CSF1R inhibition does not only increase apoptosis rate of glioma cells but also decreases their proliferation rate (26). Meanwhile, the CSF1R inhibitor-induced IGF1 secretion from macrophages through IL4 signaling could, in turn, sustains the proliferation of glioma cells (3).
d) The CSF1R inhibitor treatment does not significantly affect the number of tumor-associated macrophages, as experimental studies indicated (27). As such, the total density of tumor-associated macrophages is assumed to be constant during CSF1R inhibition treatment.
Schematic representation of how adaptation of the tumor microenvironment contributes to acquired drug resistance during CSF1R inhibition therapy of gliomas (3). Macrophages and glioma cells can interact through the CSF1/CSF1R and IGF1/IGF1R pathways. CSF1R, a critical receptor on macrophages, is a promising immunotherapeutic target in gliomas and is undergoing clinical evaluation. In the absence of CSF1R inhibition, macrophages exhibit a protumorigenic M2 phenotype that promotes glioma progression. When targeted by a CSF1R inhibitor, macrophages switch to an antitumorigenic phenotype that promotes tumor cell death. On the other hand, following prolonged treatment, elevated macrophage-derived IGF1, which is promoted by IL4 secreted by other cell types (e.g., T cells), in turn sustains the survival and growth of glioma cells.
We used a set of PDEs to model the growth, death, and migration of cancer cells and tumor-associated macrophages and the diffusion/secretion/degradation of CSF1, IGF1, and drugs. The time scale of this entire mechanism is about 200 days in accordance with the experimental study (3).
The migration of glioma cells consists of random diffusion and chemotaxis toward the chemical gradients of signals. We assumed a logistic growth rate and a linear death rate of cancer cells (28, 29). Therefore, the equation that describes the evolution of the tumor cell density (CT) is as follows:
where DT is the diffusion coefficient of glioma cells and a1 is the chemotaxis coefficient for glioma cells with respect to IGF1. 
represents the maximal carrying capacity of the cancer cells. 
and 
are the proliferation rate and the death rate, respectively, which are regulated by the IGF1 and CSF1R inhibitor. Experimental data (3, 26) have revealed that CSF1R inhibition does not only increase apoptosis of glioma cells but also decreases their proliferation. Although the CSF1R inhibitor-induced IGF1 secretion from macrophages could promote the proliferation of glioma cells (3). Therefore, the regulated proliferation rate of glioma cells is modeled by the following equation:
where 
is the basic proliferation rate of tumor cells and 
is a regulatory coefficient for the macrophage density (
) corresponding to M2 macrophages' protumorigenic function. Hill functions (30–32) for the activation and repression processes were used to describe the roles of IGF1 and CSF1R inhibitors (CSF1R_I) in regulating the proliferation rate of glioma cells as follows:
where 
and 
are Michaelis–Menten constants.
The apoptosis of glioma cells is accelerated by treatment with a CSF1R inhibitor (3, 26), as described above; therefore, the regulated death rate of glioma cells is described by the following equation:
where 
is the basic death rate of glioma cells and 
is a regulatory coefficient corresponding to M1 macrophages' antitumorigenic function. The pro-apoptotic role of the CSF1R inhibitor was modeled as
where 
is a Michaelis–Menten constant.
Tumor-associated macrophages are recruited to tumor sites by chemoattractants such as CSF1 (33, 34); this process can be modeled by the following equation:
where 
is the diffusion coefficient of the random walk of macrophages and 
is the chemotaxis coefficient of CSF1. The experimental studies have demonstrated that CSF1R inhibition did not affect the number of tumor-associated macrophages (26, 27). Therefore, the total macrophage density is assumed to be constant, as also reported in refs. (25, 35). Therefore, the net growth and death rates of macrophages are zero in the above equation. The above Eq. (4) is a simplification of the tumor-associated macrophage population which involves both M1 and M2-type macrophages. M1 macrophages exhibit antitumorigenic phenotype with enhanced phagocytosis, whereas M2 macrophages are protumorigenic and contribute to glioma growth and progression. However, CSF1R inhibition alters macrophages' phenotypes switched from M2 to M1 (27). Functions of these two phenotypes are modeled in Eqs. (2, 3), respectively. The CSF1R inhibitor induced phenotype switch of macrophages from M2 to M1 is reflected by the decrease of protumor function of M2 macrophages (
) through repression Hill function H2(CSF1R_I) and the increase of antitumor function of M1 macrophage (
) through activation Hill function H3(CSF1R_I). Furthermore, because populations of both these two types of macrophage are subject to diffusion and chemotaxis as in Eq. (4), for simplicity we do not describe the two types separately.
CSF1 is overexpressed in human glioblastomas (36) and can be synthesized and produced in both cell-surface form and secreted form (37). The secreted CSF1 by glioma cells then diffuses into the microenvironment. Therefore, the concentration of CSF1 is modeled by the following equation:
where 
represents the diffusion rate of CSF1. 
and 
are the rate of CSF1 secretion by tumor cells and the CSF1 degradation rate, respectively.
IGF1 secretion from macrophages is promoted by IL4 signaling. We model the evolution of the IGF1 concentration as follows:
where 
, 
and 
are the diffusion, secretion and degradation rates of IGF1, respectively. Following treatment with a CSF1R inhibitor, IL4 is secreted from multiple cell types, including T cells and its concentration is elevated. The experimental data (3) indicate that the numbers of many types of immune cells (e.g., neutrophils, monocytes, dendritic-like cells, B cells, helper T cells and regulatory T cells) in drug-resistant group of mice are not significantly different from that in control group. As such, we do not explicitly describe the densities of IL4-secreting cells in the model for simplicity. We assume that the rate of change in the concentration of IL4 is proportional to the concentration of the CSF1R inhibitor with a coefficient of 
; that is, the secretion of IL4 can be modeled by the following equation:
Solving this equation and substituting it into Eq. (6) yields the following equation:
where 
is the initial concentration of IL4 in the absence of treatment with the CSF1R inhibitor.
In addition, we model the diffusion and degradation of drugs (including CSF1R inhibitors and IGF1R inhibitors) as follows:
where 
is the diffusion rate of the drug and 
is the degradation rate or the uptake rate of the drug by cells.
Non-dimensionization and radial symmetry simplification
To reduce the number of parameters, we simplify the PDE model by rescaling the equations into a dimensionless form (Supplementary Text S1). Furthermore, for simplicity, we assumed a spherical model for the tumor and its microenvironment. Therefore, the domain we consider is a spheroid. Assuming radial symmetry, the cell density and cytokine concentration evolve in the domain [0, L]. The details of radial symmetry simplification are provided in Supplementary Text S2.
Initial and boundary conditions
The initial and boundary conditions are given in Supplementary Text S3. In particularly, the addition of the drug is modeled with a Dirichlet boundary condition 
.
Numerical simulation
We solve the system of the dimensionless PDEs (Eqs. Supplementary S11–S15) numerically using the pdepe solver in MATLAB R2007b (MathWorks), which automatically and dynamically selects time steps and is capable of solving PDEs with multiple temporal orders (38). The length of the simulation domain (L) is approximately 4 cm. The numerical simulation results are evaluated at locations with spatial step of 0.08 cm and time points with step of 1 day.
Parameter sensitivity analysis
Most parameters in the model are obtained from experimental measurements, whereas some others are estimated or calibrated. The details of parameter estimation were described in Supplementary Text S4. Supplementary Tables S1 and S2 list the original parameters and the dimensionless parameters, respectively, used in the simulation. To examine whether the parameter variations significantly affect the model output of interest, we implement sensitivity analyses of the most relevant parameters, i.e., AT, Aϕ, rT, dT, B1, B2, e, Dd, and 
.
A sampling-based global sensitivity analysis (39) is performed. For each parameter, we use Latin hypercube sampling to generate 1,000 samples in the range of 
, where p0 is the baseline value of the parameter. All parameters are varied simultaneously for the simulations.
For the varied parameters, we choose the killing rate and the regrowth rate of cancer cells as model outputs to evaluate the long-term effects of parameter variations on the drug response and resistance. The killing rate of cancer cells is defined as
where 
represents the average cancer cell density and t1 is the time at which the minimal average cancer cell density is achieved. Similarly, the regrowth rate of cancer cell density is defined as
where t2 is the length of time over which the simulation performed.
For each parameter, the partial correlation coefficient between the simulated model output (RK or RG) and the sampled parameter values is calculated and is defined as a global sensitivity coefficient in this study.
A local parameter sensitivity analysis (Supplementary Text S5) is also performed to investigate the dose-dependent sensitivities of parameters and to assess the combinatorial effects of the parameters on drug resistance.
Results
The evolution profiles of the tumor and microenvironment in response to CSF1R inhibitor treatment
We used our mathematical model to investigate the evolution profile of the tumor and the associated microenvironment under the condition of treatment with a CSF1R inhibitor. The drug response of a high-grade glioma tumor (with a radius of approximately 4 cm) to treatment with a CSF1R inhibitor (3 doses) was simulated (Fig. 2). As shown in Fig. 2A, at the early stage, the cancer cell density decreased dramatically after a dozen days. However, the tumor regrew after a long-term dormancy period of approximately 150 days. At day 200, the cancer cell density was restored to nearly the level that was comparable before the treatment. The changes in macrophage density and the concentrations of CSF1 and IGF1 in the tumor microenvironment during the CSF1R inhibitor treatment were also examined (Fig. 2B–D). As expected, the change in the CSF1 concentration exhibited similar pattern with the evolution of the cancer cell density, whereas the IGF1 concentration kept increasing over time. Figure 2E shows the evolution profile of the concentration of CSF1R inhibitor. In addition, Fig. 2F shows the time courses of average cancer cell density and average drug concentration. Our model observed similar regrowth of resistant gliomas under the treatment of CSF1R inhibitor as previously reported (Fig. 1B in ref. 3).
Drug response of cancer cells and adaptation of the tumor microenvironment in the context of therapeutic intervention with CSF1R inhibition. The spatio-temporal distributions of (A) cancer cell density, (B) macrophage density, (C) CSF1 concentration, (D) IGF1 concentration, and (E) drug concentration. F, Curves of time course of average cancer cell density and average drug concentration.
Key parameters affecting the dynamics of drug resistance
To determine which parameters are more sensitive and critical for the killing rate and the regrowth rate of cancer cells during treatment with a CSF1R inhibitor, we performed a global parameter sensitivity analysis (see Supplementary Text S3). Figure 3 shows the results of the global sensitivity analysis. The killing rate (Fig. 3A) and the regrowth rate (Fig. 3B) display different sensitivities with respect to variations in the parameters. The significant sensitive parameters for the killing rate are rT, dT, B2, and Dd, whereas the significant sensitive parameters for the regrowth rate are AT, rT, B2, Dd, and 
.
Parameter sensitivity analysis. A, Global sensitivity analyses of the cancer cell killing rate with respect to the model parameters. B, Global sensitivity analyses of the cancer cell regrowth rate with respect to the model parameters. Positive sensitivity indicates that increasing the corresponding parameter leads to an increase in the killing rate or the regrowth rate. (*) indicates that the sensitivity for the corresponding parameter is significantly (P < 0.05) different from zero. The drug dose was 3 in the simulation.
To confirm that the emergence of the regrowth of the cancer cells is not caused by the decreased drug concentration due to drug degradation or consumption, we examined the spatio-temporal profiles of cancer cells and their microenvironment using the model with increased drug consumption rate (10-fold of 
; Supplementary Fig. S1). The simulation results showed that increasing the drug consumption rate did not significantly influence the emergence of drug resistance, compared with Fig. 2.
We also examined the behavior of the model using different values of certain important parameters, including rT, dT, B2, and Dd (Supplementary Fig. S2). A decrease in the value of rT (growth rate of tumor cells) and an increase in dT (death rate of tumor cells), as well as a decrease in B2 (secretion rate of IGF1), resulted in the disappearance of the emergence of drug resistance. In addition, the diffusion rate of the drug (Dd) also affected the killing and regrowth rates of the resistant cancer cells.
Stochasticity in the secretion rate of IGF1 results in heterogeneous drug resistance
The secretion rate of IGF1 (B2) is a critical parameter for the macrophage-mediated drug resistance as described above (Fig. 3). Therefore, we then examined whether stochasticity in IGF1 secretion rate could result in heterogeneous drug resistance. We assumed a bimodal distribution of B2 and used the Monte Carlo method to simulate the stochastic temporal evolution of cancer cell density in a population of N glioma patients. The details of the stochastic simulation and the computation of the distribution of survival times were described in Supplementary Text S6.
Figure 4A shows the heterogeneous evolution of cancer cell density in different patients in response to drug treatment. Each curve represents one patient. The simulation demonstrated that the cancer cell density was reduced in some patients but some tumors regrew. Furthermore, Fig. 4B shows the theoretical survival curves, representing the changes of survival percentages over time, of the glioma patients with and without CSF1R inhibition treatment (red line and gray line, respectively). The predicted significant difference between these two curves and their trends are in good agreement with the preclinical data (Fig. 1D in ref. 3). Moreover, we quantified the survival percentage, 
, from the experimental data (Fig. 1D in ref. 3) and calculated the corresponding survival frequency according to 
, with results shown in Fig. 4C and D. The comparison between the simulated survival frequency and experimental data shows a good agreement between them.
Stochastic simulation and model validation with survival times. A, Stochastic time course of cancer cell density. Each curve represents one patient. B, Survival curves of patients treated with CSF1R inhibitor (solid red curve) or with no treatment (dashed gray curve). The significant difference between these two survival curves and their tendencies are in good agreement with existing preclinical data (3). C, Comparison between experimental data and prediction of the frequencies of survival times under no treatment condition. D, Comparison between experimental data and prediction of the frequencies of survival times under CSF1R inhibitor treatment condition.
Dose-dependent drug resistance
We further investigated the effect of dose on the drug response and resistance. Figure 5A displays curves for the time course of cancer cell density under conditions of treatment with different doses (0, 0.1, 0.5, 1, 2, and 3) of a CSF1R inhibitor. A low dose (e.g., less than 0.1) had a slight inhibitory effect on tumor growth. High doses killed cancer cells at the early stage but did not prevent tumor recurrence at a later stage. As the dose was increased, the rate of cancer cell death increased; however, the regrowth rate also increased due to drug resistance. A more comprehensive illustration of the “dose-temporal drug response” (Fig. 5B) demonstrated that different drug doses resulted in differential temporal evolutions of cancer cell density. The emergence of drug resistance, as well as its extent, depended on the dose of CSF1R inhibitor used.
Dose-dependent drug resistance. A, Time courses of the average cancer cell density under conditions of treatment with various doses (0, 0.1, 0.5, 1, 2, and 3) of a CSF1R inhibitor. B, Dose-temporal drug response of cancer cell density. Average density of cancer cells in the whole domain was shown.
Quantitative evaluation of the effect of combination drug therapy
To determine how the combined parameters affect the killing and regrowth rates of cancer cells during drug treatment, we performed a 2-parameter sensitivity analysis (see Supplementary Text S5). The results revealed that variations in the parameter combinations r1&d1, r1&B2, d1&B2, and d1&B2 influenced the killing rate of cancer cells more than other combinations (Supplementary Fig. S3A). For the regrowth rate, the most sensitive parameter combination was r1&B2 (Supplementary Fig. S3B). These results indicate that combination therapy using drugs that both prevent cancer cell growth and inhibit the IGF1 signaling axis might be more effective than single CSF1R inhibition therapy.
Therefore, we evaluated the results of combination therapy using a CSF1R inhibitor and an IGF1R inhibitor (Fig. 6). The effect of IGF1R inhibition was incorporated into the model using the following equation:
where K1 is a half-constant that was used to normalize the drug concentration of the IGF1R inhibitor. Therefore, K1 was set to 1 in the simulation.
Evaluation of combination therapy using a CSF1R inhibitor and an IGF1R inhibitor. A, Evolution of cancer cell densities with single or combined CSF1R and IGF1R inhibitors. B, Dose-dependent killing rate and (C) regrowth rate of cancer cells under combined CSF1R inhibition and IGF1R inhibition. D, Evaluation of the synergy between CSF1R inhibition and IGF1R inhibition using the Bliss combination index.
To quantitatively evaluate the synergy of the CSF1R-I and IGF1R-I co-treatment, a combination index based on Bliss independence (40, 41) was defined as follows:
where 
quantifies the drug efficacy of the single CSF1R inhibitor (
) or IGF1R inhibitor (
) alone at dose 
, 
in Eq. (13) is the expected efficacy of the combination drug, and 
is the actual outcome produced by the combination. Hence, the combination of CSF1R-I with IGF1R-I is synergistic if CI > 0, antagonistic if CI < 0, and additive if CI = 0.
Figure 6 demonstrates the dose-dependent effects of a combination therapy using a CSF1R inhibitor and an IGF1R inhibitor simultaneously. Figure 6A shows the time course of cancer cell density under conditions of treatment with single or combined CSF1R and IGF1R inhibitors. The combined intervention significantly reduced the resistance of the glioma to CSF1R inhibition. This prediction is consistent with experimental results (figure 7 in ref. 3) that the combination of BLZ945 (CSF1R inhibitor) and OSI906 (IGF1R inhibitor) significantly improved outcome in preclinical models. Figure 6B and C further shows the dependence of the killing and regrowth rates of the cancer cells on the drug dose. The regrowth rate decreased as the dose of IGF1R inhibitor was increased. In particular, at a dose of (1, 1), the killing rate was high, and the regrowth rate was significantly decreased. Moreover, at these doses, the combination of a CSF1R inhibitor and an IGF1R inhibitor had a synergistic effect, as indicated by a positive combination index (Fig. 6D). Furthermore, the combination of CSF1R inhibition and IGF1R inhibition exhibited a pattern of dose-dependent synergy in which increasing doses resulted in less synergy, indicating the importance of drug dosages optimization.
Discussion
The tumor microenvironment plays multiple indispensable roles in cancer progression, and as such, targeting the tumor microenvironment has been proposed as a promising strategy for treating cancers (1, 2). Macrophages are of the most abundant immune cell types in the tumor microenvironment, the accumulation of which has been associated with the progression and high grade of various tumors, including gliomas (42, 43). Preclinical experiments using genetic mouse models have been conducted to demonstrate that macrophages can be targeted by inhibiting CSF1R to regress high-grade gliomas (3). However, like other targeted therapies, acquired resistance to CSF1R inhibition also occurs following long-term treatment (3).
In this study, we developed a simplified spatio-temporal model to mechanistically describe the response of glioma cells to CSF1R inhibition and the associated dynamic adaptation of the tumor microenvironment. The model simulations of tumor evolution and drug resistance are consistent with experimental results (3). The modeling framework provides a basis for the further investigation of tumor microenvironment-mediated drug resistance, provided that patient-specific data are available in future studies. Moreover, our results revealed a dose-dependent mechanism of resistance to CSF1R inhibition. Higher doses induced a greater regrowth rate of cancer cells compared with lower doses, although the killing rate was also increased. These results imply that multiple doses should be considered in experimental designs and/or clinical trials aiming to investigate resistance mechanisms and optimize drug doses.
The global parameter sensitivity analysis varies all the parameter simultaneously and evaluates the correlation between the overall model outputs and a set of sampled parameter values. Although the local parameter sensitivity analysis evaluates the local effect of the variation of a single parameter on the model output. The local parameter sensitivity analysis could differentiate the effect of different treatment conditions. A dose-dependent parameter sensitivity was revealed by the local sensitivity analysis (Supplementary Fig. S4) when treatment with different doses of CSF1R inhibitor was modeled. At a CSF1R inhibitor dose of 1 (Supplementary Fig. S4A and S4B), the most important parameters for the killing rate were the growth rate (r1) and the death rate (d1) of the cancer cells; other parameters were much less sensitive (Supplementary Fig. S4A). As the dose was increased to 3 (Supplementary Fig. S4C and S4D), the sensitivity coefficients of the parameters decreased except for that of parameter B2, which influences the drug-induced secretion rate of IGF1 (Eq. Supplementary S14) in Supplementary Text S2). An important and interesting finding is that the relative sensitivity of the regrowth rate to the parameters r1 and d1 was altered as the drug dose increased from 1 to 3 (Supplementary Fig. S4B and S4D). Moreover, the regrowth rate was found to be more sensitive than the killing rate of cancer cells to parameter B2.
We further investigated the impact of the initial spatial distribution of macrophages on the dynamics of drug response using different initial conditions of macrophages (Supplementary Text S3). The simulation results (Supplementary Figs. S5 and S6) demonstrated that, although the total densities of macrophages in these situations were set the same, the resistance kinetics of cancer cell density and spatio-temporal distribution profiles of cells and the secreted CSF1/IGF1 were different from that in Fig. 2. In particular, the simulation suggested that the higher ratio of macrophages close to the center of the tumor region, the higher density of resistant tumor cells on 200 days. These results illustrated the importance of the spatial distribution of tumor associated macrophages in mediating the resistance of CSF1R inhibition in gliomas.
We quantitatively evaluated the effects of combination therapy using CSF1R inhibition and IGF1R inhibition. Notably, we found that this drug combination exhibited a pattern of dose-dependent synergy and that the synergy decreased as the drug doses increased. This indicates that the efficiency of high drug doses in combination was reduced due to antagonistic effects, although high doses of combined drugs might produce higher killing rates or lower regrowth rates compared with low doses of combined drugs. Considering the increased side effects induced by high doses of drugs, information on how to optimize the doses of combined drugs in combination therapy is very important, especially in the era of precision medicine. In addition, because the killing rate and the regrowth rate of cancer cells have distinct sensitivities to variations in parameters, the timing of the administration of combined drugs might be critical for their long-term treatment effect. Therefore, temporal combinations of multiple drugs with optimized timing and dosing might be more effective than the use of single drugs or simultaneous combinations of drugs (44).
The previous experimental studies have demonstrated that the operations of immune cells and other microenvironmental variables might influence the response of glioma to the treatment (45, 46). As such, besides macrophages modeled in this study, in the future study, we will consider other immune cells that are relevant in gliomas, such as dendritic cells, neutrophils, Treg cells, T effector cells, natural killer cells, etc. (14). In addition, cancer cell populations are heterogeneous and they might react differentially to immune cells. The model could be extended to include cancer stem cells to investigate the heterogeneous response of cancer cells to immunotherapies. Moreover, the in vivo data used to estimate the model parameters are indispensable to a more realistic model because the previous relevant works have showed that the slight changes in the model parameter might significantly influence the prediction results (47–49). In future work, we will develop data-driven methods to construct the signaling pathways involved in crosstalk interactions between gliomas and their microenvironments using single cell RNA-seq data and proteomics data. Furthermore, we will build multiscale agent-based models (50) of cell-microenvironment interactions to investigate the intracellular and intercellular mechanisms underlying the drug resistance of cancer cells.
In summary, in this study, we developed a PDE model to quantitatively investigate the mechanisms of interaction between glioma cells and macrophages that underlie the resistance of these tumor cells to CSF1R inhibition. Our model revealed a dose-dependent resistance mechanism and indicated that a parameter sensitivity switch is triggered by increasing the drug dose. We quantitatively evaluated the synergistic effects of the combination of CSF1R inhibition and IGF1R inhibition. This study advances our understanding of microenvironment-mediated drug resistance and has implications for the design of more effective combination drug therapies to combat drug resistance.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: X. Sun, Y. Zheng
Development of methodology: X. Sun, Y. Zheng
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Zheng
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): X. Sun, Y. Zheng
Writing, review, and/or revision of the manuscript: X. Sun, Y. Zheng, Q. Zhao, T. Zhou
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): X. Sun, Y. Zheng
Study supervision: X. Sun, Y. Zheng
Acknowledgments
X. Sun was supported by grants from the National Natural Science Foundation of China (61503419), the Guangdong Nature Science Foundation (2014A030310355, 2016A030313234), and the fund for Guangdong Provincial Key Laboratory of Orthopedics and Traumatology (2016B030301002). T. Zhou was supported by grants from the National Key Research Project of China (91530320) and the 973 Project of China (2014CB964703).
We would like to acknowledge Dr. Tao Cai at the School of Mathematical and Computational Science, Sun Yat-sen University, for valuable discussions on numerical simulation. We also thank Prof. Shicheng Su at Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University, for help with editing the manuscript.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Footnotes
Note: Supplementary data for this article are available at Molecular Cancer Therapeutics Online (http://mct.aacrjournals.org/).
- Received July 6, 2017.
- Revision received November 1, 2017.
- Accepted January 18, 2018.
- Published first February 13, 2018.
- ©2018 American Association for Cancer Research.
References
Volume 17, Issue 4
