Mathematical modeling and analysis of cancer treatment with radiation and anti-PD-L1.

In cancer treatment, radiation therapy (RT) induces direct tumor cell death due to DNA damage, but it also enhances the deaths of radiosensitive immune cells and is followed by local relapse and up-regulation of immune checkpoint ligand PD-L1. Since the binding between PD-1 and PD-L1 curtails anti-tumor immunities, combining RT and PD-L1 inhibitor, anti-PD-L1, is a potential method to improve the treatment efficacy by RT. Some experiments support this hypothesis by showing that the combination of ionizing irradiation (IR) and anti-PD-L1 improves tumor reduction comparing to the monotherapy of IR or anti-PD-L1. In this work, we create a simplified ODE model to study the order of tumor growths under treatments of IR and anti-PD-L1. Our synergy analysis indicates that both IR and anti-PD-L1 improve the tumor reduction of each other, when IR and anti-PD-L1 are given simultaneously. When giving IR and anti-PD-L1 separately, a high dosage of IR should be given first to efficiently reduce tumor load and then followed by anti-PD-L1 with strong efficacy to maintain the tumor reduction and slow down the relapse. Increasing the duration of anti-PD-L1 improves the tumor reduction, but it cannot prolong the duration that tumor relapses to the level of the control case. Under some simplification, we also prove that the model has an unstable tumor free equilibrium and a locally asymptotically stable tumor persistent equilibrium. Our bifurcation diagram reveals a transition from tumor elimination to tumor persistence, as the tumor growth rate increases. In the tumor persistent case, both anti-PD-L1 and IR can reduce tumor amount in the long term.

IL-27 in combination with anti-PD-1 can be anti-cancer or pro-cancer.

Interleukin-27 (IL-27) is known to play opposing roles in immunology. The present paper considers, specifically, the role IL-27 plays in cancer immunotherapy when combined with immune checkpoint inhibitor anti-PD-1. We first develop a mathematical model for this combination therapy, by a system of Partial Differential Equations, and show agreement with experimental results in mice injected with melanoma cells. We then proceed to simulate tumor volume with IL-27 injection at a variable dose F and anti-PD-1 at a variable dose g. We show that in some range of "small" values of g, as f increases tumor volume decreases as long as fFc(g), where Fc(g) is a monotone increasing function of g. This demonstrates that IL-27 can be both anti-cancer and pro-cancer, depending on the ranges of both anti-PD-1 and IL-27.

Different mechanisms of CD200-CD200R induce diverse outcomes in cancer treatment.

The CD200 is a cell membrane protein expressed by tumor cells, and its receptor CD200 receptor (CD200R) is expressed by immune cells including macrophages and dendritic cells. The formation of CD200-CD200R inhibits the cellular functions of the targeted immune cells, so CD200 is one type of the immune checkpoint and blockade CD200-CD200R formation is a potential cancer treatment. However, the CD200 blockade has opposite treatment outcomes in different types of cancers. For instance, the CD200R deficient mice have a higher tumor load than the wild type (WT) mice in melanoma suggesting that CD200-CD200R inhibits melanoma. On the other hand, the antibody anti-CD200 treatment in pancreatic ductal adenocarcinoma (PDAC) and head and neck squamous cell carcinoma (HNSCC) significantly reduces the tumor load indicating that CD200-CD200R promotes PDAC and HNSCC. In this work, we hypothesize that different mechanisms of CD200-CD200R in tumor microenvironment could be one of the reasons for the diverse treatment outcomes of CD200 blockade in different types of cancers. We create one Ordinary Differential Equations (ODEs) model for melanoma including the inhibition of CCL8 and regulatory T cells and the switching from M2 to M1 macrophages by CD200-CD200R to capture the tumor inhibition by CD200-CD200R. We also create another ODEs model for PDAC and HNSCC including the promotion of the polarization and suppressive activities of M2 macrophages by CD200-CD200R to generate the tumor promotion by CD200-CD200R. Furthermore, we use these two models to investigate the treatment efficacy of the combination treatment between the CD200-CD200R blockade and the other immune checkpoint inhibitor, anti-PD-1. Our result shows that different mechanisms of CD200-CD200R can induce different treatment outcomes in combination treatments, namely, only the CD200-CD200R blockade reduces tumor load in melanoma and only the anti-PD-1 and CD200 knockout decrease tumor load in PDAC and HNSCC. Moreover, in melanoma, the CD200-CD200R mainly utilizes the inhibitions on M1 macrophages and dendritic cells to inhibit tumor growth, instead of M2 macrophages.

A simple in-host model for COVID-19 with treatments: model prediction and calibration.

In this paper, we provide a simple ODEs model with a generic nonlinear incidence rate function and incorporate two treatments, blocking the virus binding and inhibiting the virus replication to investigate the impact of calibration on model predictions for the SARS-CoV-2 infection dynamics. We derive conditions of the infection eradication for the long-term dynamics using the basic reproduction number, and complement the characterization of the dynamics at short-time using the resilience and reactivity of the virus-free equilibrium are considered to inform on the average time of recovery and sensitivity to perturbations in the initial virus free stage. Then, we calibrate the treatment model to clinical datasets for viral load in mild and severe cases and immune cells in severe cases. Based on the analysis, the model calibrated to these different datasets predicts distinct scenarios: eradication with a non reactive virus-free equilibrium, eradication with a reactive virus-free equilibrium, and failure of infection eradication. Moreover, severe cases generate richer dynamics and different outcomes with the same treatment. Calibration to different datasets can lead to diverse model predictions, but combining long- and short-term dynamics indicators allows the categorization of model predictions and determination of infection severity.

G protein controls stress readiness by modulating transcriptional and metabolic homeostasis in Arabidopsis thaliana and Marchantia polymorpha.

The core G protein signaling module, which consists of Gα and extra-large Gα (XLG) subunits coupled with the Gßγ dimer, is a master regulator of various stress responses. In this study, we compared the basal and salt stress-induced transcriptomic, metabolomic and phenotypic profiles in Gα, Gß, and XLG-null mutants of two plant species, Arabidopsis thaliana and Marchantia polymorpha, and showed that G protein mediates the shift of transcriptional and metabolic homeostasis to stress readiness status. We demonstrated that such stress readiness serves as an intrinsic protection mechanism against further stressors through enhancing the phenylpropanoid pathway and abscisic acid responses. Furthermore, WRKY transcription factors were identified as key intermediates of G protein-mediated homeostatic shifts. Statistical and mathematical model comparisons between A. thaliana and M. polymorpha revealed evolutionary conservation of transcriptional and metabolic networks over land plant evolution, whereas divergence has occurred in the function of plant-specific atypical XLG subunit. Taken together, our results indicate that the shifts in transcriptional and metabolic homeostasis at least partially act as the mechanisms of G protein-coupled stress responses that are conserved between two distantly related plants.

Mathematical modeling for the combination treatment of IFN-γ and anti-PD-1 in cancer immunotherapy.

When the immune-checkpoint programmed death-1 (PD-1) binds to its ligand programmed death ligand 1 (PD-L1) to form the complex PD-1-PD-L1, this complex inactivates immune cells resulting in cell apoptosis, downregulation of immune reaction, and tumor evasion. The antibody, anti-PD-1 or anti-PD-L1, blocks the PD-1-PD-L1 complex formation to restore the functions of T cells. Combination of anti-PD-1 with other treatment shows promising in different types of cancer treatments. Interferon-gamma (IFN-γ) plays an important role in immune responses. It is mainly regarded as a pro-inflammatory cytokine that promotes the proliferation of CD8+ T cell and cytotoxic T cell, enhances the activation of Th1 cells and CD8+ T cells, and enhances tumor elimination. However, recent studies have been discovering many anti-inflammatory functions of IFN-γ, such as promotion of the PD-L1 expression, T cell apoptosis, and tumor metastasis, as well as inhibition of the immune recognition and the killing rates by T cells. In this work, we construct a mathematical model incorporating pro-inflammatory and anti-inflammatory functions of IFN-γ to capture tumor growth under anti-PD-1 treatment in the wild type and IFN-γ null mutant melanoma. Our simulation results qualitatively fit experimental data that IFN-γ null mutant with anti-PD-1 obtains the highest tumor reduction comparing to IFN-γ null mutant without anti-PD-1 and wild type tumor with anti-PD-1 therapy. Moreover, our synergy analysis indicates that, in the combination treatment, the tumor volume decreases as either the dosage of anti-PD-1 increases or the IFN-γ production efficiency decreases. Thus, the combination of anti-PD-1 and IFN-γ blockade improves the tumor reduction comparing to the monotherapy of anti-PD-1 or the monotherapy of IFN-γ blockade. We also find a threshold curve of the minimal dosage of anti-PD-1 corresponding to the IFN-γ production efficiency to ensure the tumor reduction under the presence of IFN-γ.

Mathematical Modeling and Analysis of CD200-CD200R in Cancer Treatment.

CD200 is a cell membrane protein that binds to its receptor, CD200 receptor (CD200R). The CD200 positive tumor cells inhibit the cellular functions of M1 and M2 macrophages and dendritic cells (DCs) through the CD200-CD200R complex, resulting in downregulation of Interleukin-10 and Interleukin-12 productions and affecting the activation of cytotoxic T lymphocytes. In this work, we provide two ordinary differential equation models, one complete model and one simplified model, to investigate how the binding affinities of CD200R and the populations of M1 and M2 macrophages affect the functions of the CD200-CD200R complex in tumor growth. Our simulations demonstrate that (i) the impact of the CD200-CD200R complex on tumor promotion or inhibition highly depends on the binding affinity of the CD200R on M2 macrophages and DCs to the CD200 on tumor cells, and (ii) a stronger binding affinity of the CD200R on M1 macrophages or DCs to the CD200 on tumor cells induces a higher tumor cell density in the CD200 positive tumor. Thus, the CD200 blockade would be an efficient treatment method in this case. Moreover, the simplified model shows that the binding affinity of CD200R on macrophages is the major factor to determine the treatment efficacy of CD200 blockade when the binding affinities of CD200R on M1 and M2 macrophages are significantly different to each other. On the other hand, both the binding affinity of CD200R and the population of macrophages are the major factors to determine the treatment efficacy of CD200 blockade when the binding affinities of CD200R on M1 and M2 macrophages are close to each other. We also analyze the simplified model to investigate the dynamics of the positive and trivial equilibria of the CD200 positive tumor case and the CD200 deficient tumor case. The bifurcation diagrams show that when M1 macrophages dominate the population, the tumor cell density of the CD200 positive tumor is higher than the one of CD200 deficient tumor. Moreover, the dynamics of tumor cell density change from tumor elimination to tumor persistence to oscillation, as the maximal proliferation rate of tumor cells increases.

Mathematical modeling and analysis of the heat shock protein response during thermal stress in fish and HeLa cells.

The climate change has the potential to dramatically affect species' thermal physiology and may change the ecology and evolution of species' lineages. In this work, we investigated the transition of dynamics in the heat shock response when the thermal stress approaches the upper thermal limits of species to understand how the climate change may affect the heat shock responses in ectotherms and endotherms. The heat shock protein 70, HSP70, prevents protein denaturation or misfolding under thermal stresses. When thermal stress increases, the number of misfolded proteins increases, which leads to high levels of HSP70 protein. However, when temperatures approach limits of thermal tolerance (i.e., the critical thermal maximum, CTmax, for ectotherms and the superior critical temperature, SCT, for endotherms), levels of HSP70 protein synthesis start to decrease. Thus, we hypothesized that the temperature at the first reduction of HSP abundance indicates the thermal limits of the species. In this work, we provide a mathematical model to investigate the behavior of the heat shock responses related to HSP70 protein. This model captures the dynamics of HSP70 protein and Hsp70 mRNA, in HeLa cells (i.e., representative for endotherms) and multiple species of fishes (i.e., representative for ectotherms) with different acclimation histories. Based on our hypothesis of the relationship between the HSP70 protein level and CTmax/SCT, our model provides three methods to predict the CTmax of fishes with varying acclimation temperature and the SCT of HeLa cells. The CTmax increases as the acclimation temperature increases in fishes, however the CTmax plateaus when the acclimation temperature is higher than 20°C in brook trout, a representative cool water salmonid. Our model also captures the situation that the heat shock reaction in fish is more sensitive to the heat shock temperature than HeLa cells, when the heat shock temperature is lower than the upper thermal tolerance. However, both fish and HeLa cells are sensitive to the heat shock temperature when the temperature reaches the thermal tolerance limits. Additionally, our sensitive analysis result indicates that the status of some components in the heat shock reaction changes when the temperature reaches the thermal tolerance resulting in failure in protein refolding in fish and speeding up the refolding process in HeLa cells. Mathematical analysis is also applied on a simplified mathematical model to investigate the detailed dynamics of the model, such as the steady states of the substrate, Hsp70 mRNA, and HSP70 protein, at different thermal stresses. The comparison between the original model and simplified model shows that the inhibition on HSP70 protein transcription by thermal stresses leads to the reduction in HSP70 protein at extreme temperatures.

Mathematical modeling of the eyespots in butterfly wings.

Butterfly wing color patterns are a representative model system for studying biological pattern formation, due to their two-dimensional simple structural and high inter- and intra-specific variabilities. Moreover, butterfly color patterns have demonstrated roles in mate choice, thermoregulation, and predator avoidance via disruptive coloration, attack deflection, aposematism, mimicry, and masquerade. Because of the importance of color patterns to many aspects of butterfly biology and their apparent tractability for study, color patterns have been the subjects of many attempts to model their development. Early attempts focused on generalized mechanisms of pattern formation such as reaction-diffusion, diffusion gradient, lateral inhibition, and threshold responses, without reference to any specific gene products. As candidate genes with expression patterns that resembled incipient color patterns were identified, genetic regulatory networks were proposed for color pattern formation based on gene functions inferred from other insects with wings, such as Drosophila. Particularly detailed networks incorporating the gene products, Distal-less (Dll), Engrailed (En), Hedgehog (Hh), Cubitus interruptus (Ci), Transforming growth factor-ß (TGF-ß), and Wingless (Wg), have been proposed for butterfly border ocelli (eyespots) which helps the investigation of the formation of these patterns. Thus, in this work, we develop a mathematical model including the gene products En, Hh, Ci, TGF-ß, and Wg to mimic and investigate the eyespot formation in butterflies. Our simulations show that the level of En has peaks in the inner and outer rings and the level of Ci has peaks in the inner and middle rings. The interactions among these peaks activate cells to produce white, black, and yellow pigments in the inner, middle, and outer rings, respectively, which captures the eyespot pattern of wild type Bicyclus anynana butterflies. Additionally, our simulations suggest that lack of En generates a single black spot and lack of Hh or Ci generates a single white spot, and a deficiency of TGF-ß or Wg will cause the loss of the outer yellow ring. These deficient patterns are similar to those observed in the eyespots of Vanessa atalanta, Vanessa altissima, and Chlosyne nycteis. Thus, our model also provides a hypothesis to explain the mechanism of generating the deficient patterns in these species.

Long-distance communication in Arabidopsis involving a self-activating G protein.

In plant cells, heterotrimeric G protein signaling mediates development, biotic/abiotic stress responsiveness, hormone signaling, and extracellular sugar sensing. The amount of sugar in plant cells fluctuates from nanomolar to high millimolar concentrations over time depending on changes in the light environment. Arabidopsis thaliana Regulator of G Signaling protein 1 (AtRGS1) modulates G protein activation and detects the concentration and the exposure time of sugars. This is called dose-duration reciprocity in sugar sensing and occurs through AtRGS1 internalization which is directly proportional to G protein activation. One source of sugars is from CO 2 fixation by photosynthesis. Through a simple set of experiments, we show that sugars made in cotyledons that are undergoing photomorphogenesis activate G signaling in cells distal to the nascent photosynthesis center. This occurs with sufficient speed to enable distal cells to monitor changes in photosynthetic activity in the leaves.

The kinetics in mathematical models on segmentation clock genes in zebrafish.

Somitogenesis is the process for the development of somites in vertebrate embryos. This process is timely regulated by synchronous oscillatory expression of the segmentation clock genes. Mathematical models expressed by delay equations or ODEs have been proposed to depict the kinetics of these genes in interacting cells. Through mathematical analysis, we investigate the parameter regimes for synchronous oscillations and oscillation-arrested in an ODE model and a model with transcriptional and translational delays, both with Michaelis-Menten type degradations. Comparisons between these regimes for the two models are made. The delay model has larger capacity to accommodate synchronous oscillations. Based on the analysis and numerical computations extended from the analysis, we explore how the periods and amplitudes of the oscillations vary with the degradation rates, synthesis rates, and coupling strength. For typical parameter values, the period and amplitude increase as some synthesis rate or the coupling strength increases in the ODE model. Such variational properties of oscillations depend also on the magnitudes of time delays in delay model. We also illustrate the difference between the dynamics in systems modeled with linear degradation and the ones in systems with Michaelis-Menten type reactions for the degradation. The chief concerns are the connections between the dynamics in these models and the mechanism for the segmentation clocks, and the pertinence of mathematical modeling on somitogenesis in zebrafish.

Dose-Duration Reciprocity for G protein activation: Modulation of kinase to substrate ratio alters cell signaling.

In animal cells, activation of heterotrimeric G protein signaling generally occurs when the system's cognate signal exceeds a threshold, whereas in plant cells, both the amount and the exposure time of at least one signal, D-glucose, are used toward activation. This unusual signaling property called Dose-Duration Reciprocity, first elucidated in the genetic model Arabidopsis thaliana, is achieved by a complex that is comprised of a 7-transmembrane REGULATOR OF G SIGNALING (RGS) protein (AtRGS1), a Gα subunit that binds and hydrolyzes nucleotide, a Gßγ dimer, and three WITH NO LYSINE (WNK) kinases. D-glucose is one of several signals such as salt and pathogen-derived molecular patterns that operates through this protein complex to activate G protein signaling by WNK kinase transphosphorylation of AtRGS1. Because WNK kinases compete for the same substrate, AtRGS1, we hypothesize that activation is sensitive to the AtRGS1 amount and that modulation of the AtRGS1 pool affects the response to the stimulant. Mathematical simulation revealed that the ratio of AtRGS1 to the kinase affects system sensitivity to D-glucose, and therefore illustrates how modulation of the cellular AtRGS1 level is a means to change signal-induced activation. AtRGS1 levels change under tested conditions that mimic physiological conditions therefore, we propose a previously-unknown mechanism by which plants react to changes in their environment.

A nondestructive method to estimate the chlorophyll content of Arabidopsis seedlings.

BACKGROUND: Chlorophyll content decreases in plants under stress conditions, therefore it is used commonly as an indicator of plant health. Arabidopsis thaliana offers a convenient and fast way to test physiological phenotypes of mutations and treatments. However, chlorophyll measurements with conventional solvent extraction are not applicable to Arabidopsis leaves due to their small size, especially when grown on culture dishes. RESULTS: We provide a nondestructive method for chlorophyll measurement whereby the red, green and blue (RGB) values of a color leaf image is used to estimate the chlorophyll content from Arabidopsis leaves. The method accommodates different profiles of digital cameras by incorporating the ColorChecker chart to make the digital negative profiles, to adjust the white balance, and to calibrate the exposure rate differences caused by the environment so that this method is applicable in any environment. We chose an exponential function model to estimate chlorophyll content from the RGB values, and fitted the model parameters with physical measurements of chlorophyll contents. As proof of utility, this method was used to estimate chlorophyll content of G protein mutants grown on different sugar to nitrogen ratios. CONCLUSION: This method is a simple, fast, inexpensive, and nondestructive estimation of chlorophyll content of Arabidopsis seedlings. This method lead to the discovery that G proteins are important in sensing the C/N balance to control chlorophyll content in Arabidopsis.

A shadow detector for photosynthesis efficiency.

Plants tolerate large variations in the intensity of the light environment by controlling the efficiency of solar to chemical energy conversion. To do this, plants have a mechanism to detect the intensity, duration, and change in light as they experience moving shadows, flickering light, and cloud cover. Sugars are the primary products of CO2 fixation, a metabolic pathway that is rate limited by this solar energy conversion. We propose that sugar is a signal encoding information about the intensity, duration and change in the light environment. We previously showed that the Arabidopsis heterotrimeric G protein complex including its receptor-like Regulator of G signaling protein, AtRGS1, detects both the concentration and the exposure time of sugars (Fu et al., 2014. Cell 156: 1084-1095). This unique property, designated dose-duration reciprocity, is a behavior that emerges from the system architecture / system motif. Here, we show that another property of the signaling system is to detect large changes in light while at the same time, filtering types of fluctuation in light that do not affect photosynthesis efficiency. When AtRGS1 is genetically ablated, photosynthesis efficiency is reduced in a changing- but not a constant-light environment. Mathematical modeling revealed that information about changes in the light environment is encoded in the amount of free AtRGS1 that becomes compartmentalized following stimulation. We propose that this property determines when to adjust photosynthetic efficiency in an environment where light intensity changes abruptly caused by moving shadows on top of a background of light changing gradually from sun rise to sun set and fluctuating light such as that caused by fluttering leaves.

The role of the cytokines IL-27 and IL-35 in cancer.

The cancer-immune interaction is a fast growing field of research in biology, where the goal is to harness the immune system to fight cancer more effectively. In the present paper we review recent work of the interaction between T cells and cancer. CD8+ T cells are activated by IL-27 cytokine and they kill tumor cells. Regulatory T cells produce IL-35 which promotes cancer cells by enhancing angiogenesis, and inhibit CD8+ T cells via TGF-ß production. Hence injections of IL-27 and anti-IL-35 are both potentially anti-tumor drugs. The models presented here are based on experimental mouse experiments, and their simulations agree with these experiments. The models are used to suggest effective schedules for drug treatment.

Mathematical modeling of Interleukin-35 promoting tumor growth and angiogenesis.

Interleukin-35 (IL-35), a cytokine from the Interleukin-12 cytokine family, has been considered as an anti-inflammatory cytokine which promotes tumor progression and tumor immune evasion. It has also been demonstrated that IL-35 is secreted by regulatory T cells. Recent mouse experiments have shown that IL-35 produced by cancer cells promotes tumor growth via enhancing myeloid cell accumulation and angiogenesis, and reducing the infiltration of activated CD8[Formula: see text] T cells into tumor microenvironment. In the present paper we develop a mathematical model based on these experimental results. We include in the model an anti-IL-35 drug as treatment. The extended model (with drug) is used to design protocols of anti-IL-35 injections for treatment of cancer. We find that with a fixed total amount of drug, continuous injection has better efficacy than intermittent injections in reducing the tumor load while the treatment is ongoing. We also find that the percentage of tumor reduction under anti-IL-35 treatment improves when the production of IL-35 by cancer is increased.

Mathematical modeling of interleukin-27 induction of anti-tumor T cells response.

Interleukin-12 is a pro-inflammatory cytokine which promotes Th1 and cytotoxic T lymphocyte activities, such as Interferon-[Formula: see text] secretion. For this reason Interleukin-12 could be a powerful therapeutic agent for cancer treatment. However, Interleukin-12 is also excessively toxic. Interleukin-27 is an immunoregulatory cytokine from the Interleukin-12 family, but it is not as toxic as Interleukin-12. In recent years, Interleukin-27 has been considered as a potential anti-tumor agent. Recent experiments in vitro and in vivo have shown that cancer cells transfected with IL-27 activate CD8+ T cells to promote the secretion of anti-tumor cytokines Interleukin-10, although, at the same time, IL-27 inhibits the secretion of Interferon-[Formula: see text] by CD8+ T cells. In the present paper we develop a mathematical model based on these experimental results. The model involves a dynamic network which includes tumor cells, CD8+ T cells and cytokines Interleukin-27, Interleukin-10 and Interferon-[Formula: see text]. Simulations of the model show how Interleukin-27 promotes CD8+ T cells to secrete Interleukin-10 to inhibit tumor growth. On the other hand Interleukin-27 inhibits the secretion of Interferon-[Formula: see text] by CD8+ T cells which somewhat diminishes the inhibition of tumor growth. Our numerical results are in qualitative agreement with experimental data. We use the model to design protocols of IL-27 injections for the treatment of cancer and find that, for some special types of cancer, with a fixed total amount of drug, within a certain range, continuous injection has better efficacy than intermittent injections in reducing the tumor load while the treatment is ongoing, although the decrease in tumor load is only temporary.

The effect of time delay in a two-patch model with random dispersal.

We consider a two-patch model for a single species with dispersal and time delay. For some explicit range of dispersal rates, we show that there exists a critical value τc for the time delay τ such that the unique positive equilibrium of the system is locally asymptotically stable for τ ∈[0,τc) and unstable for τ > τc .

The role of CD200-CD200R in tumor immune evasion.

CD200 is a cell membrane protein that interacts with CD200 receptor (CD200R) of myeloid lineage cells. During tumor initiation and progression, CD200-positive tumor cells can interact with M1 and M2 macrophages through CD200-CD200R-compex, and downregulate IL-10 and IL-12 productions secreted primarily by M2 and M1 macrophages, respectively. In the tumor microenvironment, IL-10 inhibits the activation of cytotoxic T lymphocytes (CTL), while IL-12 enhances CTL activation. In this paper, we used a system approach to determine the combined effect of CD200-CD200R interaction on tumor proliferation by developing a mathematical model. We demonstrate that blocking CD200 on tumor cells may have opposite effects on tumor proliferation depending on the "affinity" of the macrophages to form the CD200-CD200R-complex with tumor cells. Our results help understanding the complexities of tumor microenvironment.