Modeling Macrophage Polarization and Its Effect on Cancer Treatment Success

Show more

1. Introduction

1.1. Immune Cell Polarization

The differentiation and polarization of certain immune cells into pro-inflammatory and anti-inflammatory cells confer the immune system with versatility to exhibit a strong, but controlled, response against invading pathogens and foreign antigens. This is made possible by the initial generation of a strong inflammatory response that is subsequently regulated and attenuated by an anti-inflammatory response once the invading agents have been destroyed. A phenotypic switch by immune cells during an infection and after its resolution makes it possible for the immune system to regulate its own activity and return to a state of homeostasis [1] .

When healthy cells mutate and become cancerous, immune cells such as natural killer cells and cytotoxic T lymphocytes (CTL) try to eliminate these anomalous cells. They do so by infiltrating the tumor site and releasing proteins that destroy the cancer cell membrane (Perforin), and release enzymes that lead to cancer cell apoptosis (Granzyme B) [2] . Other immune cells, such as macrophages, are capable of directly phagocytosing the cancer cells. Macrophages that attempt to kill tumor cells via contact-dependent mechanisms or through molecule secretion are classified as M1 macrophages. M2 macrophages tend to exhibit a pro-tumor phenotype characterized by the release of immunosuppressive and angiogenic agents that allow a tumor to survive and grow. It is believed that the M1-to-M2 macrophage concentration ratio can tip the balance in favor of tumor destruction if this ratio is high enough, or in favor of tumor survival if this ratio is close to zero [3] . Macrophage repolarization rates influence the M1-to-M2 ratio which, in turn, could be used to predict tumor size [4] . Similarly, naive helper T cells differentiate into anti-tumor Th1 helper cells or into pro-tumor Th2 helper cells depending on the relative concentrations of anti-tumor and pro-tumor cytokines found in the tumor microenvironment [5] . Consequently, the Th1-to-Th2 helper cell ratio also affects the likelihood of tumor destruction or survival [6] . A high Th1-to-Th2 ratio increases the likelihood of tumor destruction, whereas a ratio close to zero increases the likelihood of tumor survival [7] .

Experimental work has shown that it is possible to reprogram M2 macrophages to develop an M1 phenotype by increasing the environmental concentration of anti-tumor cytokines. Monocytes and naïve helper T cells develop an M1 and Th1 phenotype, respectively, if the concentration of anti-tumor cytokinesINF-γ and TNF-α is high. Given this phenotypic plasticity of immune cells, it is no surprise that cancer cells have evolved ways to hijack the mechanisms of cell polarization for their own advantage. Cancer cells can escape immune destruction by releasing pro-tumor cytokines, such as TGF-β, that decrease the M1-to-M2 macrophage ratio and the Th1-to-Th2 helper cell ratio. What makes this hijacking relevant, in the context of cancer treatments, is that macrophages and helper T cells are found in the tumor site in relatively high concentrations [8] [9] [10] [11] . The high proportion of M2 macrophages and Th2 helper cells in the tumor microenvironment makes them an important target of cancer therapies [12] [13] [14] [15] [16] .

1.2. Positive Feedback Loops Perpetuate Cell Polarization

A tumor is a complex dynamical system, and its survival depends on a diverse set of signaling networks characterized by cytokine-driven positive feedback loops that can reinforce the anti-tumor phenotype or the pro-tumor phenotype of tumor-infiltrating immune cells. For example, M1 macrophages secrete IL-12 which leads to the differentiation of immature helper T cells into Th1 cells. Th1 cells secrete IFN-ϒ which reinforces the M1 macrophage phenotype. This positive feedback loop perpetuates the M1 and Th1 anti-tumor polarization of these cells, which can lead to tumor destruction. On the other hand, M2 macrophages secrete IL-4 and IL-6 [17] , [18] which lead to the differentiation of immature helper T cells into Th2 cells. Th2 cells secrete IL-4 [19] which reinforces the M2 macrophage phenotype. This positive feedback loop perpetuates the M2 and Th2 pro-tumor polarization, leading to tumor escape.

More complex immune cell interactions exist. M2 macrophages and Th2 cells secrete TGF-β which converts naïve helper T cells into pro-tumor regulatory T cells (Tregs) [20] and B cells into pro-tumor regulatory B cells (Bregs). Tregs and Bregs secrete IL-10 and TGF-β, which reinforces the Th2 phenotype. Th2 cells then reinforce the M2 phenotype and, as a result, the pro-tumor activity of all these cells is perpetuated. Moreover, by releasing IL-10 and TGF-β, Tregs and Bregs reinforce each other’s pro-tumor characteristics. The presence of positive feedback loops can lead to switch-like dynamical behavior and bi-stability [21] , [22] . Self-perpetuating positive feedback loops and switch-like bi-stability can make a tumor robust with respect to external perturbations, such as a cancer treatment. Positive feedback loops generate attracting dynamical states from which it is difficult to escape without a strong external perturbation. It is for this reason that combination treatments have become a promising approach to treat cancer. Figure 1 illustrates the complex interactions between cancer cells and the immune system, between the immune cells, and the positive feedback loops that characterize such interactions. A thorough review of these interactions can be found in [23] .

This article is organized as follows. In Section 2 we describe the mathematical model that was formulated to investigate potential ways to break up pro-tumor feedback loops that lead to treatment failure. We also used the model to assess the relative effectiveness of various combination treatments. In Section 3 we describe the predictions of the model, including the combination treatments that were found to be robust with respect to the level of immune response strength, the tumor size at the start of treatment, treatment resistance, and immune cell polarization. In Section 4 we discuss the implications of the model predictions and comment on ways that the model can be improved. We conclude by highlighting the importance of increasing the M1-to-M2 ratio and the Th1-to-Th2

Figure 1. Illustration of the tumor-immune system interactions.

ratio to boost the anti-tumor immune response, in conjunction with a reduction of TGF-β-driven angiogenesis. This approach is predicted to lead to treatment synergy and robustness, and to the effective disruption of the pro-tumor cytokine-driven signaling networks that lead to tumor survival.

2. Mathematical Model

2.1. Model Description

In [24] , Fernandez and Soto-Ortiz expanded an ordinary differential equation (ODE) model of drug resistance [25] to include tumor angiogenesis stimulated by TGF-β and immunosuppression exerted by regulatory T cells. Our model expands [24] to include cytokine-driven feedback loops that lead to the polarization of macrophages and helper T cells into anti-tumor cells or pro-tumor cells. The model includes ODEs describing the change in the concentration of M1 and M2 macrophages, Th1 and Th2 helper cells, and the concentration of anti-tumor cytokines IL-1β, IL-12, TNF-α and IFN-γ and pro-tumor cytokines IL-4, IL-10, IFN-α and TGF-β. We modeled IL-2 as an anti-tumor and pro-tumor cytokine because IL-2 can stimulate NK cells [26] , CTL [27] and Tregs [28] .

The ODE equations were coded in Scilab (http://www.scilab.org/) and were solved by using a built-in 4th-order explicit Runge-Kutta method with fixed step size. This numerical method has a fast rate of convergence of O(h^{4}) and guarantees a stable computation time. Many of the model parameter values were obtained from published literature pertaining to experimental or mathematical modeling work involving various murine and human cancers, as well as infection models that describe a pathogen-immune system interaction that is similar to a tumor-immune system interaction.

It is known that the phenotypes of macrophages vary widely and, thus, they cannot be classified as being purely anti-tumor or purely pro-tumor [29] . Macrophages tend to exhibit a heterogeneity of phenotypes that sometimes overlap. To simplify the development and analysis of the model, we assumed that all the differentiated macrophages have an anti-tumor or pro-tumor phenotype, and that monocytes can differentiate into M1 or M2 macrophages. We also assumed that when an anti-tumor or pro-tumor macrophage repolarizes by changing its phenotype, it immediately adopts the opposite phenotype without delay and without passing through an intermediate state.

The Th1 and Th2 phenotypes of helper T cells tend to be terminally-differentiated states [23] . Repolarization between the Th1 and Th2 states seldom occurs under natural conditions. Experimentally, Th1 and Th2 cells can be repolarized under specific experimental conditions [30] . Hence, our model assumes that once naive helper T cells differentiate and polarize into a Th1 or Th2 state, they maintain that terminal phenotype throughout their existence. The Th1-to-Th2 cell ratio cannot be modulated directly through repolarization of Th1 or Th2 cells. However, the Th1-to-Th2 cell ratio can be modulated indirectly through macrophage repolarization which then polarizes naïve helper T cells, or by direct infusion of polarized helper T cells.

2.2. Model Variables and Initial Conditions

The equations of the model are presented in Appendix A. The ODEs that describe the population dynamics of treatment-sensitive and treatment-resistant cancer cells include terms that represent cancer cell destruction by NK cells, CTL, M1 macrophages and Th1 helper T cells. Chemotherapy kills cancer cells and immune system cells. We did not incorporate an anti-tumor humoral response. In the cytokine equations, we included terms that represent antagonistic interactions between certain cytokines. For example, an increase in the concentration of IL-4 and IL-10 decreases the production rate of IL-12 by M1 macrophages and by Th1 helper cells. Similarly, an increase in the concentration of IL-12 decreases the production rate of IL-4 by M2 macrophages and by Th2 helper cells. Antagonistic effects between IFN-α and TNF-α, and between IFN-α and IFN-γ were also included. Supplementary Table S1 in Appendix B lists the definitions of the model variables and their units.

To investigate the joint effect of immunosuppression, tumor angiogenesis and immune cell polarization on treatment success, we considered three different scenarios at the start of treatment. These scenarios are listed below in order from most favorable to least favorable:

1) Low immunosuppression, low angiogenesis and no pro-tumor polarization.

2) High immunosuppression, high angiogenesis and no pro-tumor polarization.

3) High immunosuppression, high angiogenesis and high pro-tumor polarization.

To simulate scenario 1, the initial conditions were chosen by assuming a low concentration of Tregs, TGF-β and activated endothelial cells, and by assuming that no M2 macrophages and no Th2 cells are initially present. To simulate scenario 2, the initial conditions were chosen by assuming a high concentration of Tregs, TGF-β and activated endothelial cells, but no M2 macrophages and no Th2 cells are initially present. To simulate scenario 3, the initial conditions were chosen to reflect a high concentration of Tregs, TGF-β, activated endothelial cells, M2 macrophages and Th2 cells. Supplementary Table S2 in Appendix B lists the initial conditions that represent each of these scenarios. Many of these initial conditions are the same that were used in [24] . For the case of an initially high pro-tumor cell polarization, we used the M2 macrophage and Th2 cell concentrations predicted in the tumor microenvironment according to the mathematical model presented in [4] , and the IL-4 and IL-10 concentrations used in [31] .

2.3. Simulated Treatments

Supplementary Table S3 in Appendix B lists all the monotherapies that we considered, including a description of the dose, frequency and the time required for administration. These monotherapies served as the building blocks of the combination treatments that we simulated. ODEs were used to simulate the treatments as intravenous injections of constant infusion rate, or as capsules taken orally. The infusion rate terms that are listed were obtained by following the procedures described in [32] and we refer to them as the regular rates. Since M2 macrophages and Tregs are found in the tumor site at high concentrations, the rates of infusion of M1 macrophages and of Th1 helper cells were assigned the same value as the rate of infusion of NK cells. These M1 and Th1 infusion rates make it possible to significantly increase the M1-to-M2 macrophage ratio and the Th1-to-Th2 helper cell ratio. Doing so allowed us to quantify the effect of modulating these polarization ratios on treatment success.

As was done in [24] , we assumed that wild-type cancer cells become resistant to Irinotecan chemotherapy in a dose-dependent manner. We also considered a hypothetical chemotherapy drug to which no cancer cells become resistant. Moreover, we assumed that in all simulations there were 35 KRAS-mutant cancer cells which are resistant to the monoclonal antibodies Panitumumab and Cetuximab. An objective of our project was to identify treatment combinations that are robust with respect to treatment resistance and that can eliminate the wild-type and the resistant cancer cells.

2.4. Model Parameters and Treatment Robustness

To investigate the effect of cytokine concentration on immune cell polarization and tumor growth, we used a partial differential equation (PDE) model of a granuloma that develops in response to a lung infection [31] . That model describes the immune reaction to an infection of Mycobacterium tuberculosis. It includes PDEs that simulate the spatiotemporal dynamics of pro- and anti-tumor alveolar macrophages, helper T cells, CTL and dendritic cells. The infection model also includes PDEs that describe the dynamics of the cytokines that are responsible for the polarization of alveolar macrophages and naïve helper T cells. The immune reaction against this bacterial infection is very similar to an immune reaction against cancer cells: the innate, adaptive and humoral immune responses are activated. Like the lung bacteria-immune system interaction, the same anti-tumor and pro-tumor cytokines play a role in the polarization of tumor-associated macrophages and helper T cells. A series of signaling pathways are turned on and off over time that lead to an initial inflammatory response characterized by a high concentration of M1 macrophages, Th1 cells and anti-tumor cytokines followed by an anti-inflammatory response characterized by a high concentration of M2 macrophages, Th2 cells and pro-tumor cytokines. Once the bacterial infection is resolved, the immune system returns to a state of homeostasis. In the case of cancer development, the immune response is unable to eliminate the tumor and becomes anergic. A high concentration of pro-tumor cells and cytokines remains in the tumor site, making the microenvironment highly immunosuppressive and the eradication of the tumor more difficult.

Supplementary Table S4 in Appendix B lists the parameter values that were used in the simulations. We incorporated into the ODE model [24] the cytokine-driven immune cell polarization and the antagonistic interactions of various cytokines simulated in the PDE infection model [31] . However, we assumed a homogeneous concentration of immune cells and cytokines in the tumor microenvironment. Hence, we used only the ordinary derivative terms of the PDE infection model. We also assumed that the rates of immune cell growth, polarization and cytokine production by polarized macrophages and helper T cells during a Mycobacterium tuberculosis infection are similar to those in a tumor-immune system interaction. The parameter values that we estimated were chosen based on similarities between two processes, or between cytokines that exert similar functions, as explained in [24] . The parameter values in [31] were converted to appropriate units for use in our model. For example, the units of concentration of the cytokines were converted from mg/L to IU/L by referring to the published specific activity of each cytokine. See [33] for details on how to perform these conversions. Cytokines whose specific activity was not available in the published literature were assigned an average specific activity of 1.3 × 10^{10} IU/g. Moreover, the density of macrophages, helper T cells and cytokines were reported in [31] in g/cm^{3}. We converted the density of cells to cell/L and the density of cytokines to IU/L. The wet weight of an epithelial tumor having a volume of 1 cm^{3} is approximately 1 gram and contains approximately 10^{8} cancer cells [34] . Although different types of cells vary in size, we assumed a uniform sizefor cancer cells, macrophages and helper T cells and that there are 10^{8} cells per gram.

The parameters d, l and s together determine the level of strength of a patient’s immune response D by cytotoxic T lymphocytes (CTL), as defined by Equation (28). The 3 levels of immune response strength that we considered were the same as those defined in [32] . To better understand the effect of cytokine-driven pro-tumor cell polarization on treatment success, we considered nine cases. They are listed below in order from best-case scenario to worst-case scenario at the start of treatment:

1) Strong Immune Response: $\left(d=2.1,l=1.1,s=5\times {10}^{-3}\right)$

a) Best case: Low immunosuppression, low angiogenesis and no polarization.

b) High immunosuppression, high angiogenesis and no polarization.

c) High immunosuppression, high angiogenesis and high pro-tumor polarization.

2) Moderate Immune Response: $\left(d=1.6,l=1.4,s=8\times {10}^{-3}\right)$

a) Low immunosuppression, low angiogenesis and no polarization.

b) High immunosuppression, high angiogenesis and no polarization.

c) High immunosuppression, high angiogenesis and high pro-tumor polarization.

3) Weak Immune Response: $\left(d=1.3,l=2,s=4\times {10}^{-2}\right)$

a) Low immunosuppression, low angiogenesis and no polarization.

b) High immunosuppression, high angiogenesis and no polarization.

c) Worst case: High immunosuppression, high angiogenesis and high pro-tumor polarization.

To ensure the safety in the clinic of all the simulated treatments, we followed the safety criteria described in [24] . In all the simulations, a complete response (CR) to a cancer treatment was defined as a tumor size at the end of the treatment period that is less than or equal to the diffusion-limited value of 1 × 10^{6} tumor cells. A partial response (PR) to treatment is described by a tumor that remains larger than 1 × 10^{6} cells, but that by the end of the treatment period is smaller than at the start of treatment. The no response (NR) classification applies when tumor size remains the same, or if the tumor becomes larger than it was at the start of treatment, by the time the treatment period ends. A cancer treatment that leads to a complete response at the three levels of immune response strength, for tumor sizes of up to 10^{9} cells, and that kills cancer cells that are resistant to Irinotecan and Panitumumab, was defined as being robust with respect to these perturbations. We defined a treatment that eliminated a tumor despite very low initial M1-to-M2 and Th1-to-Th2 ratios as being robust with respect to pro-tumor cell polarization.

3. Results

3.1. Strong Immune Response

We first considered the case of initially low immunosuppression, low angiogenesis and no initial pro-tumor cell polarization (refer to Supplementary Table S2 in Appendix B for the initial conditions that represent this scenario) coupled with a strong immune response. This scenario is not the typical one in the case of tumor growth. Tumors evade detection and destruction by reducing their antigenicity and by suppressing the immune system. Consequently, the immune response is attenuated and fails to eliminate the tumor. However, when we assumed the favorable conditions above, this led to the destruction of the tumor without the need for treatment, as can be seen in Figure 2. At steady-state, the M1-to-M2 macrophage polarization ratio is no more than one order of magnitude greater than 1, while the Th1-to-Th2 helper cell polarization ratio approaches 1. This result represents immune cell homeostasis after the cancer cells have been eliminated. The model predictions for this first scenario highlight the importance of reducing immunosuppression and tumor angiogenesis to allow the immune system to naturally eliminate the tumor.

Figure 3 illustrates the case of initially high immunosuppression, high angiogenesis and no initial pro-tumor immune cell polarization. The immune system is unable to eliminate the tumor without treatment despite its ability to produce a strong CTL anti-tumor response. The M1-to-M2 and Th1-to-Th2 polarization ratios both remain close to zero due to the high initial concentration of M2 macrophages and of Th2 helper cells and throughout the simulation time. This occurs due to self-reinforcing positive feedback loops that lead to concentrations of pro-tumor cytokines that remain several orders of magnitude higher than the concentrations of the anti-tumor cytokines. The lopsided imbalance of these cytokine concentration perpetuates animmune cell polarization that favors tumor survival.

Under the initial conditions described above, the model predicts that a5-cycle Sunitinib + NK cell treatment given at their regular rates will fail to reduce the size of the tumor despite the significant reduction of the Treg population caused by the Sunitinib injections. This treatment fails due to the high concentration of immunosuppressive TGF-β at the start of the treatment and that remains high throughout the treatment period. Additionally, the model predicts that a Fresolimumab monotherapy consisting of 20 injections at the regular rate will also fail to eliminate the tumor. Similarly, 20 injections of M1 macrophages administered concurrently with 20 injections of Th1 cells at their regular rates is not sufficient to reduce the size of the tumor. An important result of this particular simulation is that injecting M1 macrophages and Th1 helper cells concurrently at their regular rates may not be sufficient to significantly alter the M1-to-M2 ratio and the Th1-to-Th2 ratio, if immunosuppression and angiogenesis are initially high. Other measures must be taken to modify these ratios, such as increasing the treatment dosage or disrupting additional signaling pathways. For example, one M1macrophage injection + one Fresolimumab (anti-TGF-β) injection administered concurrently at their regular rates eliminated the tumor in approximately 90 days. This success was due in part to the strong CTL response that the immune system is capable of exhibiting in this scenario, and the complementary

Figure 2. Without treatment, elimination of a tumor under a strong CTL response occurs whenever immunosuppression, tumor angiogenesis and immune cell polarization are low. The M1-to-M2 ratio and the Th1-to-Th2 ratio both increase, leading to an increased production of the anti-tumor cytokines IL-12, TNF-α and IFN-γ by M1 macrophages and Th1 helper cells.

Figure 3. Initially high immunosuppression, high angiogenesis and no cell polarization under a weak immune response leads to a low M1-to-M2 ratio and a low Th1-to-Th2 ratio without treatment. An M1 macrophage + Fresolimumab combination treatment eliminates the tumor in approximately 90 days by boosting the M1-to-M2 and Th1-to-Th2 ratios.

nature at work of the two treatment modalities: a boosting of the CTL response coupled with a reduction of tumor angiogenesis and immunosuppression exerted by TGF-β.

Figure 4 illustrates some unfavorable and favorable outcomes when assuming an initially high immunosuppression, high angiogenesis and high pro-tumor immune cell polarization. The combined administration of Irinotecan + Panitumumab eliminates some types of cancer cells, but not all. The wild-type cancer cells are killed by Irinotecan, by Panitumumab and by the immune system. The Irinotecan-resistant cells are eliminated by Panitumumab and by the immune system. However, over 50 Irinotecan injections and 25 Panitumumab injections must be administered, making this treatment impractical. Moreover, the Panitumumab-resistant cancer cells (KRAS-mutant) will survive. The result of this combination treatment will be a tumor consisting of Panitumumab-resistant cancer cells.

A more practical combination treatment that eliminates the tumor consists of giving 9 M1 macrophage injections + 5 Fresolimumab injections both administered at 2 times their regular rate. The need for an increased dosage was due to the high pro-tumor immune cell polarization that existed at the start of treatment. The model predicts that a monotherapy consisting of17 M1 macrophage injections administered at three times the regular rate will also eliminate the tumor. These predictions highlight the importance of reducing tumor angiogenesis and immunosuppression as a treatment strategy. These results also suggest that a lack of reduction of tumor angiogenesis can be compensated by significantly increasing the M1-to-M2 macrophage polarization ratio by increasing the infusion rates of M1 macrophages into the tumor site. We surmise that a high M1-to-M2 polarization ratio leads to a rate of cancer cell destruction by immune cells that is greater than the rate of cancer cell replication, even when this replication rate is enhanced by angiogenesis.

3.2. Moderate Immune Response

In the no treatment case with low immunosuppression, low angiogenesis and no initial immune cell polarization favorable to the tumor, the immune system initially decreases the size of the tumor. However, due to its moderate response, the immune system is unable to eliminate all the cancer cells and the tumor grows again to its maximum carrying capacity. With treatment, there were multiple treatment combinations that eliminated the tumor. A moderate immune response, without an initial pro-tumor immune cell polarization, made it possible for a single injection of Irinotecan or a single injection of Panitumumab administered at their regular rate to eliminate the tumor, including the Irinotecan-resistant and Panitumumab-resistant cells. This result is shown in Figure 5. An Irinotecan monotherapy can kill the wild-type and KRAS-mutant cells, while the Irinotecan-resistant cells are eliminated by the immune system. The moderate immune response with initially low immunosuppression, low angiogenesis and without initial pro-tumor polarization made it possible to eliminate the

Figure 4. The KRAS-mutant cancer cells survive Panitimumab treatment but are eliminated by an M1 macrophage + Fresolimumab combination therapy in 172 days.

Figure 5. Tumor escape without treatment under a moderate response. Irinotecan monotherapy and Panitumumab monotherapy eliminate the wild-type and treatment-resistant cancer cells.

tumor by giving a single M1 macrophage or Th1 cell injection at a lower infusion rate (1 × 10^{7} cells L^{−}^{1} day^{−1}) than the regular rate (5.627 × 10^{8} cells L^{−1} day^{−1}). However, a monotherapy consisting of 5 Sunitinib cycles administered at the regular rate had no effect on tumor growth. A reasonable explanation for this outcome is the fact that in this scenario, the Treg population was already low at the start of treatment and, hence, the Treg-targeted Sunitinib injections served no useful purpose.

In the case of high immunosuppression, high angiogenesis and no initial pro-tumor immune cell polarization, there were several combinations treatments that led to a complete response, as can be seen in Figure 6:

・ 6 M1 macrophage injections + 3 Fresolimumab injections given at their regular rates

・ 8 Th1 cell injections + 4 Fresolimumab injections given at their regular rates

・ 1 cycle of Sunitinib administered at its regular ratein combination with 3 injections of Fresolimumab given at 1.28 times its regular rate

・ 9 Fresolimumab injections administered at 1.45 times the regular rate

Compared to Fresolimumab monotherapy, the combination treatments eliminated the tumor in a significantly shorter amount of time.

Figure 7 shows the case of high immunosuppression, high angiogenesis and high pro-tumor immune cell polarization at the start of treatment. A Sunitinib + Fresolimumab treatment failed to eliminate the tumor even when administered at5 times their regular rates. A successful treatment consisted of administering14M1 macrophage injections together with 7 Fresolumumab injections, each given at two times their regular rate.Administering18M1 macrophage injections given at 3 times the regular rate also eliminated the tumor. The need for an increased dosage is due to the negative synergistic effect of high immunosuppression, angiogenesis and pro-tumor cell polarization at the start of treatment.

3.3. Weak Immune Response

In the case of low immunosuppression, low angiogenesis and no initial pro-tumor cell polarization (see Figure 8), an Irinotecan treatment at its regular rate eliminated the tumor with four injections. The chemotherapy-resistant cancer cells and the KRAS mutants were also eliminated due to a tumor microenvironment that was non-immunosuppressive, non-angiogenic and that contained no pro-tumor cells at the start of the treatment. As a result, the anti-tumor cytokines eventually increase in concentration, leading to M1-to-M2 and Th1-to-Th2 ratios that are representative of immune system homeostasis. In contrast, a treatment consisting of 22 Panitumumab injections, given at the regular rate, reduces the number of wild-type cancer cells, but the number of KRAS-mutant cells increases to the maximum carrying capacity, leading to treatment failure. Other successful treatments that eliminated all cancer cells consisted of a single injection of either Fresolimumab, M1 macrophages or Th1 helper cells given at their regular rate.

Figure 6. Combination treatments lead to a significantly shorter time to elimination of the tumor compared to Fresolimumab monotherapy.

Figure 7. Pro-tumor cell polarization at the start of treatment requires increased doses of M1 Macrophages and Fresolimumab to eliminate a tumor.

Figure 8. Irinotecan chemotherapy can eliminate the wild-type and mutant cancer cells if there are no pro-tumor polarized immune cells at the start of treatment, despite a weak CTL response. The cytokine plot shows that the concentration of pro-tumor cytokines IL-4, IL-10 and IFN-α released by M2 macrophages and Th2 helper cells decrease over time, leading to the success of Irinotecan chemotherapy.

In the case of high immunosuppression, high angiogenesis and no initial pro-tumor immune cell polarization, there were several combinations treatments that led to tumor elimination. The most notable case, shown in Figure 9, is the combination of 8 M1 macrophage injections coupled with 4 Fresolimumab injections, both given at their regular rates.

The M1 macrophage + Fresolimumab combination treatment was one of several treatments that were identified as being robust with respect to the level of immune response strength by CTL, with respect to the initial tumor size and with respect to resistance to Irinotecan and Panitumumab. These robust treatments, which are administered at their regular rates, are listed in Table 1.

In the worst-case scenario of high immunosuppression, high angiogenesis and high pro-tumor cell polarization with a weak immune response, there were no treatment combinations that could be administered at feasible rates that led to tumor elimination. Therefore, none of the treatments listed in Table 1 were found to be robust with respect to pro-tumor cell polarization. The failure of these promising treatments shows the negative impact on treatment success of the various positive feedback loops that reinforce pro-tumor immune cell polarization.

The lack of treatment robustness with respect to cell polarization, when assuming a weak immune response, motivated the authors to investigate the effect of directly disrupting the cytokine-driven feedback loops that are responsible for treatment failure. In general, when human cancers are diagnosed, they are well established and have already developed strong immunosuppressive mechanisms [35] . Therefore, we focused exclusively on the worst-case scenario of a weak immune response, high immunosuppression, angiogenesis, and cell polarization at the start of treatment. We noted that Fresolimumab reduces the concentration of free TGF-β in the tumor site. Therefore, Fresolimumab has the potential to reduce immunosuppression, tumor angiogenesis, and to some extent, to repolarize M2 macrophages back to the M1 phenotype. However, IL-10 and IL-4 are also responsible for pro-tumor cell polarization. Administering anti-IL-10 monotherapy biweekly or in combination with cell-based treatments did not lead to a decrease in tumor size because IL-4 plays a role that is similar to that of IL-10. Therefore, we simulated a gene knockout experiment that reduces by 99% the production of IL-4 and IL-10 by M2 macrophages and Th2 cells. We achieved this by decreasing the values of the parameters ${\lambda}_{I4M2}$ , ${\lambda}_{I4T2}$ , ${\lambda}_{I10M2}$ and ${\lambda}_{I10T2}$ by 99%. By doing so, the concentrations of IL-4 and IL-10 remained very low throughout the entire simulation time.

The gene knockout simulation showed that the tumor is eliminated when 4 M1 macrophage injections are administered together with 2Fresolimumab injections at their regular injection rates (see Figure 10). The tumor was also eliminated when 7 Th1 helper cell injections were administered together with 4Fresolimumab injections at their regular injection rates. These results show that neutralizing IL-4 and IL-10 should be part of an M1 macrophage + Fresolimumab

Figure 9. An M1 macrophage + Fresolimumab combination is an example of a treatment that is robust with respect to immune response strength, initial tumor size and treatment resistance.

Figure 10. Reducing IL-4 and IL-10 production by 99% via a gene-knockout experiment makes an M1 macrophage + Fresolimumab treatment, as well as a Th1 helper cell + Fresolimumab treatment, robust with respect to pro-tumor cell polarization.

or a Th1 cell + Fresolimumab treatment. This multi-pronged strategy confers these combination treatments with a robustness with respect to anti-tumor cell polarization, making them the most effective protocols that were identified.

3.4. Sensitivity Analysis

We proceeded to investigate the extent to which the predicted treatment outcomes depend on the values of the model parameters. To that end, we conducted a comprehensive sensitivity analysis of the weak immune response case with

Table 1. Treatments robust to immune strength level, tumor size and resistance.

high immunosuppression, high angiogenesis and no cell polarization. The analysis consisted of decreasing (and increasing) each parameter value by 5% and keeping track of the percent change in the predicted tumor size at steady state. Table 2 lists some of the parameters that were varied and the resulting percent change in the tumor size. The results of this analysis indicated that the model is sensitive to changes to three parameters: p_{1} representing the maximum rate of TGF-β production by hypoxic tumor cells, b_{K} representing the proliferation rate of angiogenic endothelial cells, and K_{max} which is the rate at which TGF-β stimulates tumor growth. The sensitivity analysis indicates that the predictions of the model depend on the rate of tumor angiogenesis and growth. Therefore, we ensured that the pharmacokinetic and pharmacodynamical properties of TGF-β are simulated appropriately.

Of note is the fact that our model is significantly less sensitive to the production rate of IL-4 and IL-10 by M2 macrophages and by Th2 cells. This means that the predicted treatment outcomes are not affected by slight changes to the rate of production and decay of these cytokines.

4. Discussion

In this work, we focused on increasing the M1-to-M2 macrophage and the Th1-to-Th2 helper cell ratios by injecting anti-tumor polarized cells (M1 macrophages and Th1 helper cells) and by disrupting pro-tumor positive feedback loops driven by TGF-β, IL-4 and IL-10. An alternative approach to increase the M1-to-M2 macrophage and the Th1-to-Th2 helper cell ratios is to inject anti-tumor cytokines, such as IL-12 and IFN-γ, to reinforce the anti-tumor positive feedback loops that polarize macrophages and naïve helper T cells into the M1 and Th1 phenotypes, respectively. We did not include in our model an immunosuppressive positive feedback loop that involves the humoral immune system.

Table 2. Sensitivity of the final tumor size to a 5% change to the parameter values.

It is known that Tregs and Bregs can reinforce each other’s pro-tumor phenotype through secretion of TGF-β and IL-10. We plan to include additional feedback loops into the model that involve the humoral immune response to assess their effect on treatment robustness. We will also consider additional sources of IL-10, such as Tregs and Bregs [36] .

The cell-cell, cell-cytokine and cytokine-cytokine interactions were modeled according to previously published models by using first-order kinetics and by using Hill functions to account for rate saturation. Although there is not a prescribed way to simulate a given interaction, it is important to simulate an interaction in the most biologically-realistic manner. Therefore, it is worth checking whether the model predictions are sensitive to the approach taken when modeling certain interactions. In the future, we plan to undertake such an analysis by introducing Hill Functions of different orders to quantify their effect on treatment success.

In our model, we chose to make the maximum tumor size and tumor angiogenesis depend on TGF-β due to its multiple pro-tumor functions [37] . However, there are other cytokines that drive angiogenesis and that are also immunosuppressive, including the Vascular Endothelial Growth Factor (VEGF) [38] [39] . We will investigate whether the treatments listed in Table S4 remain robust when we let VEGF play the role of TGF-β. We will simulate injecting the monoclonal antibody bevacizumab (Avastin) to neutralize VEGF and will compare the effectiveness of an anti-VEGF treatment versus an anti-TGF-β treatment.

In the weak immune response case with high immunosuppression, high angiogenesis, and high pro-tumor cell polarization, it is possible to eliminate the tumor by administering a smaller number of Th1 helper cell injections. However, the drawback is that it would take longer to eliminate the tumor (over 300 days). This means that we would need to keep IL-4 and IL-10 at minimal concentration during this entire treatment period. This may not be advisable, especially if there are serious safety concerns of a prolonged reduced concentration of IL-4 and IL-10. For example, it has been shown that IL-10-deficient mice develop lethal uncontrolled inflammation of the intestine [40] . IL-10 stimulates the cytotoxicity of CTL [41] and a deficiency of IL-10 can lead to spontaneous tumor development [42] . Due to safety concerns, we plan to simulate concurrent injections of anti-IL-4 and anti-IL-10 instead of reducing the parameter values that represent IL-4 and IL-10 production by M2 macrophages and by Th2 cells. By doing so, we will be able to simulate a temporary and safer neutralization of IL-4 and IL-10, and will assess the effect of this cytokine blockade on treatment success.

Treatment effectiveness does not necessarily transfer from one cancer type to another. For example, it is quite possible that a treatment that can eliminate a slow-growing tumor will fail to stop the growth of a more aggressive tumor, such as glioblastoma. This aspect will be considered when we parametrize our model to identify the treatments that are most likely to eradicate a specific type of cancer. The purpose of our modeling project was primarily to assess the extent to which cytokine-driven feedback loops that perpetuate macrophage and helper T cell polarization into a pro-tumor phenotype determine treatment outcome. Second, we wanted to quantify the extent to which disrupting such feedback loops increases the likelihood of treatment success. We conducted a relative comparison between treatments of the time required for tumor elimination, the number of required injections and the required dose. The model results predicted that injecting M1 macrophages or Th1 cells, administering anti-TGF-β to reduce angiogenesis and simultaneously reducing the production of IL-4 and IL-10 is a promising strategy to eliminate a tumor.

5. Conclusion

Cytokine-driven feedback loops play an essential role in determining the steady-state dynamics of tumor-immune cell interactions, and the likelihood of treatment success. M1 macrophage + Fresolimumab and Th1 cell + Fresolimumab combination treatments were found to be robust with respect to the level of immune response strength, the initial tumor size and treatment resistance. The model predicts that treatments that simultaneously decrease tumor angiogenesis and boost the concentration of M1 macrophages or Th1 cells will be most effective, since they boost the M1-to-M2 and Th1-to-Th2 ratios. These treatments can be made robust with respect to pro-tumor immune cell polarization if coupled with antibodies that neutralize IL-4 and IL-10. Novel cancer treatments based on IL-10 and IL-4 antibodies could pave the way for tumor elimination despite a worst-case scenario at the start of treatment consisting of a highly immunosuppressive, angiogenic and polarized tumor microenvironment.

Acknowledgements

The authors would like to thank the BUILD-PODER program at California State University-Northridge for providing the computational equipment and financial support that made this work possible.

Authors’ Contributions

LS conceived and managed the project. LS designed the computational experiments. LS and VM conducted the simulations. LS and VM analyzed the results. LS and VM contributed to manuscript preparation.

Conflicts of Interest

The authors declare that there are no conflicts of interest associated with the publication of this article.

Appendix A

Model Equations

Equations (1)-(27) represent the ODEs of the mathematical model that we used to investigate the role of cytokine-driven polarization of immune cells on treatment success. Equation (28) was defined in [32] to set the level of strength of the anti-tumor immune response by cytotoxic T lymphocytes (CTL). It was expanded further by [24] and [25] to consider KRAS mutant and Irinotecan-resistant tumor cells and the immunosuppressive effect of TGF-β on the killing rate of tumor cells by CTL. We have expanded this term further to include the immunosuppressive effect of IL-10 on CTL. Equation (29) defines the cytokine-dependent polarization of monocytes and macrophages into an M1 phenotype and of naïve helper T cells into a Th1 phenotype. Similarly, Equation 30 defines the cytokine-dependent polarization of monocytes and macrophages into an M2 phenotype and of naïve helper T cells into a Th2 phenotype.

Wild-type Cancer Cells (sensitive to all the treatments):

$\begin{array}{c}\frac{\text{d}{T}_{w}}{\text{d}t}={a}_{w}{T}_{w}\left(1-\frac{{T}_{w}+{T}_{cp}+{T}_{i}}{{T}_{K}+{g}_{2}E}\right)-\frac{\mu {C}_{1}}{{K}_{M}+{C}_{1}}{T}_{w}-\left(c+\xi \frac{A}{A+{h}_{1}}\right)\left(\text{e}{}^{-{\lambda}_{T}R}\right)N{T}_{w}-D{T}_{w}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left({K}_{T}+{K}_{AT}A\right)\left(\frac{{T}_{w}}{{T}_{w}+{\alpha}_{1}{T}_{cp}}\right)\left(1-\text{e}{}^{-{\delta}_{T}\left({C}_{1}+{C}_{2}\right)}\right){T}_{w}-\psi A{T}_{w}-{\delta}_{M1}{M}_{1}{T}_{w}-{\delta}_{T1}{T}_{1}{T}_{w}\end{array}$ (1)

KRAS-Mutant Cancer Cells (resistant to Panitumumab and Cetuximab):

$\begin{array}{c}\frac{\text{d}{T}_{cp}}{\text{d}t}={a}_{m}{T}_{cp}\left(1-\frac{{T}_{w}+{T}_{cp}+{T}_{i}}{{T}_{K}+{g}_{2}E}\right)-\left(c+\xi \frac{A}{A+{h}_{1}}\right)\left(\text{e}{}^{-{\lambda}_{T}R}\right)N{T}_{cp}-D{T}_{cp}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{T}\left(\frac{{T}_{w}}{{T}_{w}+{\alpha}_{1}{T}_{cp}}\right)\left(1-\text{e}{}^{-{\delta}_{TR}\left({C}_{1}+{C}_{2}\right)}\right){T}_{cp}-{\delta}_{M1}{M}_{1}{T}_{cp}-{\delta}_{T1}{T}_{1}{T}_{cp}\end{array}$ (2)

Chemotherapy-Resistant Cancer Cells (resistant to Irinotecan):

$\begin{array}{c}\frac{\text{d}{T}_{i}}{\text{d}t}={a}_{m}{T}_{i}\left(1-\frac{{T}_{w}+{T}_{cp}+{T}_{i}}{{T}_{K}+{g}_{2}E}\right)+\frac{\mu {C}_{1}}{{K}_{M}+{C}_{1}}{T}_{w}-\left(c+\xi \frac{A}{A+{h}_{1}}\right)\left(\text{e}{}^{-{\lambda}_{T}R}\right)N{T}_{i}-D{T}_{i}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left({K}_{T}+{K}_{AT}A\right)\left(\frac{{T}_{w}}{{T}_{w}+{\alpha}_{1}{T}_{cp}}\right)\left(1-\text{e}{}^{-{\delta}_{TR}{C}_{2}}\right){T}_{i}-\psi A{T}_{i}-{\delta}_{M1}{M}_{1}{T}_{i}-{\delta}_{T1}{T}_{1}{T}_{i}\end{array}$ (3)

Natural Killer Cells:

$\begin{array}{c}\frac{\text{d}N}{\text{d}t}=eC-fN-\left(p+{p}_{A}\frac{A}{A+{h}_{1}}\right)N\left({T}_{w}+{T}_{cp}+{T}_{i}\right)+\frac{{p}_{N}N{I}_{2}}{{g}_{N}+{I}_{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{N}\left(1-{\text{e}}^{-{\delta}_{N}\left({C}_{1}+{C}_{2}\right)}\right)N+{v}_{NK}\end{array}$ (4)

Cytotoxic T Lymphocytes:

$\begin{array}{c}\frac{\text{d}L}{\text{d}t}=-\frac{\theta mL}{\theta +{I}_{2}}+j\frac{{T}_{w}+{T}_{cp}+{T}_{i}}{k+{T}_{w}+{T}_{cp}+{T}_{i}}L-qL\left({T}_{w}+{T}_{cp}+{T}_{i}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({r}_{1}N+{r}_{2}C\right)\left({T}_{w}+{T}_{cp}+{T}_{i}\right)-\frac{u{L}^{2}R{I}_{2}}{\kappa +{I}_{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{L}\left(1-{\text{e}}^{-{\delta}_{L}\left({C}_{1}+{C}_{2}\right)}\right)L+\frac{{p}_{{I}_{2}}L{I}_{2}}{{g}_{{I}_{2}}+{I}_{2}}+{v}_{CTL}\end{array}$ (5)

Circulating T Lymphocytes in the Blood:

$\frac{\text{d}C}{\text{d}t}=\alpha -\beta C-{K}_{C}\left(1-{\text{e}}^{-{\delta}_{C}\left({C}_{1}+{C}_{2}\right)}\right)C$ (6)

Activated Endothelial Cells:

$\frac{\text{d}E}{\text{d}t}={b}_{K}E\left(1-\frac{E}{\frac{{K}_{max}B}{{g}_{1}+B}}\right)-{d}_{K}E{\left({T}_{w}+{T}_{cp}+{T}_{i}\right)}^{\frac{2}{3}}$ (7)

Regulatory T Cells

$\begin{array}{c}\frac{\text{d}R}{\text{d}t}=wC-{u}_{R}R+\left(\frac{{p}_{R}R{I}_{2}}{{g}_{R}+{I}_{2}}\right)\left(\frac{B}{B+{g}_{TGF\beta}}\right)+{\lambda}_{Treg}{T}_{0}\left(\frac{B}{B+{K}_{TGF\beta}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{h}_{R}\left(1-{\text{e}}^{-{\lambda}_{R}S}\right)R-{K}_{R}\left(1-{\text{e}}^{-{\delta}_{R}\left({C}_{1}+{C}_{2}\right)}\right)R\end{array}$ (8)

Chemotherapy Drug #1 (Irinotecan―Dose-dependent resistance emerges)

$\frac{\text{d}{C}_{1}}{\text{d}t}=-{\gamma}_{1}{C}_{1}+{v}_{C1}$ (9)

Chemotherapy Drug #2 (Hypothetical―No resistance emerges)

$\frac{\text{d}{C}_{2}}{\text{d}t}=-{\gamma}_{2}{C}_{2}+{v}_{C2}$ (10)

Interleukin-2:

$\frac{\text{d}{I}_{2}}{\text{d}t}=-{\mu}_{{I}_{2}}{I}_{2}+\varphi C+\frac{\omega L{I}_{2}}{\varsigma +{I}_{2}}+{\lambda}_{I2T1}{T}_{1}+{v}_{I2}$ (11)

Panitumumab and Cetuximab (monoclonal antibodies):

$\frac{\text{d}A}{\text{d}t}=-\eta A-\lambda \left({T}_{w}+{T}_{cp}+{T}_{i}\right)\left(\frac{A}{A+{h}_{2}}\right)+{v}_{A}$ (12)

TGF-β:

$\frac{\text{d}B}{\text{d}t}={p}_{1}\frac{{\left({T}_{w}+{T}_{cp}+{T}_{i}\right)}^{2}}{{b}_{1}^{2}+{\left({T}_{w}+{T}_{cp}+{T}_{i}\right)}^{2}}-{u}_{1}B-{b}_{2}FB$ (13)

Sunitinib (receptor tyrosine kinase inhibitor):

$\frac{\text{d}S}{\text{d}t}=-{\eta}_{S}S+{v}_{S}$ (14)

Fresolimumab (anti-TGF-β):

$\frac{\text{d}F}{\text{d}t}=-{\eta}_{F}F-{b}_{3}BF+{v}_{F}$ (15)

M1 Macrophages (anti-tumor):

$\begin{array}{c}\frac{\text{d}{M}_{1}}{\text{d}t}={\lambda}_{M0}\frac{{\epsilon}_{1}}{{\epsilon}_{1}+{\epsilon}_{2}}{M}_{0}+{\lambda}_{M1}\frac{{\epsilon}_{1}}{{\epsilon}_{1}+{\epsilon}_{2}}{M}_{2}-{\lambda}_{M2}\frac{{\epsilon}_{2}}{{\epsilon}_{2}+{\epsilon}_{1}}{M}_{1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{C}\left(1-{\text{e}}^{-{\delta}_{C}\left({C}_{1}+{C}_{2}\right)}\right){M}_{1}-{d}_{M1}{M}_{1}+{v}_{M1}\end{array}$ (16)

M2 Macrophages (pro-tumor):

$\begin{array}{c}\frac{\text{d}{M}_{2}}{\text{d}t}={\lambda}_{M0}\frac{{\epsilon}_{2}}{{\epsilon}_{2}+{\epsilon}_{1}}{M}_{0}-{\lambda}_{M1}\frac{{\epsilon}_{1}}{{\epsilon}_{1}+{\epsilon}_{2}}{M}_{2}+{\lambda}_{M2}\frac{{\epsilon}_{2}}{{\epsilon}_{2}+{\epsilon}_{1}}{M}_{1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{C}\left(1-{\text{e}}^{-{\delta}_{C}\left({C}_{1}+{C}_{2}\right)}\right){M}_{2}-{d}_{M2}{M}_{2}\end{array}$ (17)

Th1 Helper T Cells (anti-tumor):

$\begin{array}{c}\frac{\text{d}{T}_{1}}{\text{d}t}={\lambda}_{T1M1}{T}_{0}\frac{{M}_{1}}{{M}_{1}+{K}_{M1}}\frac{{I}_{12}}{{I}_{12}+{K}_{I12}}\frac{{K}_{I10}}{{K}_{I10}+{I}_{10}}+{\lambda}_{TI2}\frac{{I}_{2}}{{I}_{2}+{K}_{I2}}{T}_{1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{K}_{C}\left(1-{\text{e}}^{-{\delta}_{C}\left({C}_{1}+{C}_{2}\right)}\right){T}_{1}-{d}_{T1}{T}_{1}+{v}_{T1}\end{array}$ (18)

Th2 Helper T Cells (pro-tumor):

$\frac{\text{d}{T}_{2}}{\text{d}t}={\lambda}_{T2}{T}_{0}\frac{{M}_{2}}{{M}_{2}+{K}_{M2}}\frac{{I}_{4}}{{I}_{4}+{K}_{I4}}\frac{{K}_{T1}}{{K}_{T1}+{T}_{1}}-{K}_{C}\left(1-{\text{e}}^{-{\delta}_{C}\left({C}_{1}+{C}_{2}\right)}\right){T}_{2}-{d}_{T2}{T}_{2}$ (19)

Interleukin-1β (anti-tumor):

$\frac{\text{d}{I}_{1\beta}}{\text{d}t}=\left({\lambda}_{1\beta M1}{M}_{1}+{\lambda}_{I1\beta T1}{T}_{1}\right)\left(\frac{{K}_{I\alpha}}{{K}_{I\alpha}+{I}_{\alpha}}\right)-{d}_{I1\beta}{I}_{1\beta}$ (20)

Interleukin-4 (pro-tumor):

$\frac{\text{d}{I}_{4}}{\text{d}t}=\left({\lambda}_{I4M2}{M}_{2}+{\lambda}_{I4T2}{T}_{2}\right)\left(\frac{{K}_{I12}}{{K}_{I12}+{I}_{12}}\right)-{d}_{I4}{I}_{4}$ (21)

Interleukin-10 (pro-tumor):

$\frac{\text{d}{I}_{10}}{\text{d}t}={\lambda}_{I10M2}{M}_{2}+{\lambda}_{I10T2}{T}_{2}-{b}_{4}{A}_{I10}{I}_{10}-{d}_{I10}{I}_{10}$ (22)

Interleukin-12 (anti-tumor):

$\frac{\text{d}{I}_{12}}{\text{d}t}=\frac{{\lambda}_{I12M1}{M}_{1}+{\lambda}_{I12T1}{T}_{1}}{\left(1+\frac{{I}_{4}}{{K}_{I4}}\right)\left(\frac{1+{I}_{10}}{{K}_{I10}}\right)}-{d}_{I12}{I}_{12}+{v}_{I12}$ (23)

Tumor Necrotic Factor-α (anti-tumor):

$\frac{\text{d}{T}_{\alpha}}{\text{d}t}=\frac{{\lambda}_{T\alpha M1}{M}_{1}+{\lambda}_{T\alpha T1}{T}_{1}}{\left(1+\frac{{I}_{10}}{{K}_{I10}}\right)\left(\frac{1+{I}_{\alpha}}{{K}_{I\alpha}}\right)}-{d}_{T\alpha}{T}_{\alpha}$ (24)

Interferon-α (pro-tumor):

$\frac{\text{d}{I}_{\alpha}}{\text{d}t}={\lambda}_{I\alpha {M}_{2}}{M}_{2}\left(\frac{{K}_{I1\beta}}{{K}_{I1\beta}+{I}_{1\beta}}\right)-{d}_{I\alpha}{I}_{\alpha}$ (25)

Interferon-γ (anti-tumor):

$\frac{\text{d}{I}_{\gamma}}{\text{d}t}={\lambda}_{I\gamma T1}{T}_{1}\frac{{K}_{I\alpha}}{{K}_{I\alpha}+{I}_{\alpha}}-{d}_{I\gamma}{I}_{\gamma}+{v}_{I\gamma}$ (26)

Anti-Interleukin-10 (anti-tumor):

$\frac{\text{d}{A}_{I10}}{\text{d}t}=-{\eta}_{AI10}{A}_{I10}-{b}_{5}{I}_{10}{A}_{I10}+{v}_{AI10}$ (27)

Immune Response Strength by CTL:

$D=\frac{d}{\left(1+\frac{B}{{S}_{1}}\right)\left(1+\frac{{I}_{10}}{{T}_{I10}}\right)}\cdot \frac{{\left(\frac{L}{{T}_{w}+{T}_{cp}+{T}_{i}}\right)}^{l}}{s+{\left(\frac{L}{{T}_{w}+{T}_{cp}+{T}_{i}}\right)}^{l}}$ (28)

Polarization into Anti-Tumor Cells Driven by IFN-γ and TNF-α:

${\epsilon}_{1}=\left({\lambda}_{MI\gamma}\frac{{I}_{\gamma}}{{I}_{\gamma}+{K}_{I\gamma}}+{\lambda}_{MT\alpha}\frac{{T}_{\alpha}}{{T}_{\alpha}+{K}_{T\alpha}}\right)\left(\frac{{K}_{I10}}{{K}_{I10}+{I}_{10}}\right)\left(\frac{{K}_{TGF\beta}}{{K}_{TGF\beta}+B}\right)$ (29)

Polarization into Pro-Tumor Cells Driven by IL-4 and IL-10:

${\epsilon}_{2}=\left({\lambda}_{MI4}\frac{{I}_{4}}{{I}_{4}+{K}_{I4}}+{\lambda}_{MI10}\frac{{I}_{10}}{{I}_{10}+{K}_{I10}}\right)\left(\frac{B}{B+{K}_{TGF\beta}}\right)$ (30)

Appendix B

The Model Variables, Initial Conditions and Monotherapies

Table S1. The variables of the model.

Table S2. The initial conditions of the model.

Table S3. Monotherapy dose and frequency.

Table S4. The parameters of the model.

References

[1] Andersen, M.H. (2015) Immune Regulation by Self-Recognition: Novel Possibilities for Anticancer Immunotherapy. Journal of the National Cancer Institute, 107, djv154-djv154.

https://doi.org/10.1093/jnci/djv154

[2] Nouroz, F., Bibi, F., Noreen, S. and Masood, N. (2016) Natural Killer Cells Enhance the Immune Surveillance of Cancer. The Egyptian Journal of Medical Human Genetics, 17, 149-154.

https://doi.org/10.1016/j.ejmhg.2015.08.006

[3] Zhang, M.m et al. (2014) A high M1/M2 Ratio of Tumor-Associated Macrophages Is Associated with Extended Survival in Ovarian Cancer Patients. Journal of Ovarian Research, 7, 19.

https://doi.org/10.1186/1757-2215-7-19

[4] Den Breems, N.Y. and Eftimie, R. (2016) The Re-Polarisation of M2 and M1 Macrophages And Its Role on Cancer Outcomes. Journal of Theoretical Biology, 390, 23-39.

https://doi.org/10.1016/j.jtbi.2015.10.034

[5] Bergmann, C., Van Hemmen, J.L. and Segel, L.A. (2001) Th1 or Th2: How an Appropriate T Helper Response Can Be Made. Bulletin of Mathematical Biology, 63, 405-431.

https://doi.org/10.1006/bulm.2000.0215

[6] Kogan, Y., Agur, Z. and Elishmereni, M. (2013) A Mathematical Model for the Immunotherapeutic Control of the Th1/Th2 Imbalance in Melanoma. Discret. Contin. Dyn. Syst. - Ser. B, 18, 1017-1030.

[7] Shurin, M.R., Lu, L., Kalinski, P., Stewart-Akers, A.M. and Lotze, M.T. (1999) Th1/Th2 Balance in Cancer, Transplantation and Pregnancy. Springer Seminars In immunopathology, 21, 339-359.

https://doi.org/10.1007/BF00812261

[8] Wang, X.L., et al. (2017) High Infiltration of CD68-Tumor Associated Macrophages, Predict Poor Prognosis in Kazakh Esophageal Cancer Patients. International Journal of Experimental Pathology, 10, 10282-10292.

[9] Yafei, Z., Jun, G. and Guolan, G. (2016) Correlation between Macrophage Infiltration and Prognosis of Ovarian Cancer—A Preliminary Study. Biomedical Research, 27, No. 2.

[10] Tanaka, A. and Sakaguchi, S. (2017) Regulatory T Cells in Cancer Immunotherapy. Cell Research, 27, 109-118.

https://doi.org/10.1038/cr.2016.151

[11] Tosolini, M., et al. (2011) Clinical Impact of Different Classes of Infiltrating T Cytotoxic and Helper Cells (Th1, th2, treg, th17) in Patients with Colorectal Cancer. Cancer Research, 71, 1263-1271.

https://doi.org/10.1158/0008-5472.CAN-10-2907

[12] Zheng, X., et al. (2017) Redirecting Tumor-Associated Macrophages to Become Tumoricidal Effectors as a Novel Strategy for Cancer Therapy. Oncotarget, 8, 48436-48452.

https://doi.org/10.18632/oncotarget.17061

[13] Genard, G., Lucas, S. and Michiels, C. (2017) Reprogramming of Tumor-Associated Macrophages with Anticancer Therapies: Radiotherapy versus Chemo- and Immunotherapies. Frontiers in Immunology, 8, 828.

https://doi.org/10.3389/fimmu.2017.00828

[14] Fridlender, Z.G., et al. (2013) Using Macrophage Activation to Augment Immunotherapy of Established Tumours. British Journal of Cancer, 108, 1288-1297.

https://doi.org/10.1038/bjc.2013.93

[15] Brown, J.M., Recht, L. and Strober, S. (2017) The Promise of Targeting Macrophages in Cancer Therapy. Clinical Cancer Research, 23, 3241-3250.

https://doi.org/10.1158/1078-0432.CCR-16-3122

[16] Noy, R. and Pollard, J.W. (2014) Tumor-Associated Macrophages: From Mechanisms to Therapy. Immunity, 41, 49-61.

https://doi.org/10.1016/j.immuni.2014.06.010

[17] La Flamme, A.C., Kharkrang, M., Stone, S., Mirmoeini, S., Chuluundorj, D. and Kyle, R. (2012) Type II-Activated Murine Macrophages Produce IL-4. PLoS One, 7, e46989.

https://doi.org/10.1371/journal.pone.0046989

[18] Neurath, M.F., Fichtner-Feigl, S. and Kesselring, R. (2015) Intestinal Inflammation and Cancer of the Gastrointestinal Tract. In: Mucosal Immunology, Elsevier, Amsterdam, 1761-1775.

https://doi.org/10.1016/B978-0-12-415847-4.00091-4

[19] Romagnani, S. (1992) Type 1 T Helper and Type 2 T Helper Cells: Functions, Regulation and Role in Protection and Disease. International Journal of Clinical and Laboratory Research, 21, 152-158.

https://doi.org/10.1007/BF02591635

[20] Chen, W., et al. (2003) Conversion of Peripheral CD4+CD25- Naive T Cells to CD4+CD25+ Regulatory T Cells by TGF-Beta Induction of Transcription Factor Foxp3. Journal of Experimental Medicine, 198, 1875-1886.

https://doi.org/10.1084/jem.20030152

[21] Mitrophanov, A.Y. and Groisman, E.A. (2008) Positive Feedback in Cellular Control Systems. Bioessays, 30, 542-555.

https://doi.org/10.1002/bies.20769

[22] Ingolia, N.T. and Murray, A.W. (2007) Positive-Feedback Loops as a Flexible Biological Module. Current Biology, 17, 668-677.

https://doi.org/10.1016/j.cub.2007.03.016

[23] Burkholder, B., et al. (2014) Tumor-Induced Perturbations of Cytokines and Immune Cell Networks. Biochimica et Biophysica Acta (BBA)—Reviews on Cancer, 1845, 182-201.

https://doi.org/10.1016/j.bbcan.2014.01.004

[24] Fernandez, M., Zhou, M. and Soto-Ortiz, L. (2018) A Computational Assessment of the Robustness of Cancer Treatments with Respect to Immune Response Strength, Tumor Size and Resistance. International Journal of Tumor Therapy, 7, 1-19.

[25] Sameen, S., Barbuti, R., Milazzo, P., Cerone, A., Del Re, M. and Danesi, R. (2016) Mathematical Modeling of Drug Resistance Due to KRAS Mutation in Colorectal Cancer. Journal of Theoretical Biology, 389, 263-273.

https://doi.org/10.1016/j.jtbi.2015.10.019

[26] Jounaidi, Y., Cotten, J.F., Miller, K.W. and Forman, S.A. (2017) Tethering IL2 to Its Receptor IL2Rβ Enhances Antitumor Activity and Expansion of Natural Killer NK92 Cells. Cancer Research, 77, 5938-5951.

https://doi.org/10.1158/0008-5472.CAN-17-1007

[27] Kim, H., Kwon, B. and Sin, J.-I. (2013) Combined Stimulation of IL-2 and 4-1BB Receptors Augments the Antitumor Activity of E7 DNA Vaccines by Increasing Ag-Specific CTL Responses. PLoS One, 8, e83765.

https://doi.org/10.1371/journal.pone.0083765

[28] Liu, C., Workman, C.J. and Vignali, D.A.A. (2016) Targeting Regulatory T Cells in Tumors. The FEBS Journal, 283, 2731-2748.

https://doi.org/10.1111/febs.13656

[29] Martinez, F.O. and Gordon, S. (2014) The M1 and M2 Paradigm of Macrophage Activation: Time for Reassessment. F1000 Prime Reports, 6, 13.

https://doi.org/10.12703/P6-13

[30] Magombedze, G., Reddy, P.B.J., Eda, S. and Ganusov, V.V. (2013) Cellular and Population Plasticity of Helper CD4+ T Cell Responses. Frontiers in Physiology, 4, 206.

https://doi.org/10.3389/fphys.2013.00206

[31] Hao, W., Schlesinger, L.S. and Friedman, A. (2016) Modeling Granulomas in Response to Infection in the Lung. PLoS One, 11, e0148738.

https://doi.org/10.1371/journal.pone.0148738

[32] de Pillis, L.G., Savage, H. and Radunskaya, A.E. (2014) Mathematical Model of Colorectal Cancer with Monoclonal Antibody Treatments. British Journal of Medicine and Medical Research, 4, No. 16.

[33] Kim, K.S., Cho, G. and Jung, I.H. (2014) Optimal Treatment Strategy for a Tumor Model under Immune Suppression. Computational and Mathematical Methods in Medicine, 2014, 206287.

https://doi.org/10.1155/2014/206287

[34] Del Monte, U. (2009) Does the Cell Number 10(9) Still Really Fit one Gram of Tumor Tissue? Cell Cycle, 8, 505-506.

https://doi.org/10.4161/cc.8.3.7608

[35] Whiteside, T.L. (2010) Inhibiting the Inhibitors: Evaluating Agents Targeting Cancer Immunosuppression. Expert Opinion on Biological Therapy, 10, 1019-1035.

https://doi.org/10.1517/14712598.2010.482207

[36] Facciabene, A., Motz, G.T. and Coukos, G. (2012) T-Regulatory Cells: Key Players in Tumor Immune Escape and Angiogenesis. Cancer Research, 72, 2162-2171.

https://doi.org/10.1158/0008-5472.CAN-11-3687

[37] Neuzillet, C., et al. (2015) Targeting the TGFβ Pathway for Cancer Therapy. Pharmacology & Therapeutics, 147, 22-31.

https://doi.org/10.1016/j.pharmthera.2014.11.001

[38] Niu, G. and Chen, X. (2010) Vascular Endothelial Growth Factor as an Anti-Angiogenic Target for Cancer Therapy. Current Drug Targets, 11, 1000-1017.

https://doi.org/10.2174/138945010791591395

[39] Carmeliet, P. (2005) VEGF as a Key Mediator of Angiogenesis in Cancer. Oncology, 69, 4-10.

https://doi.org/10.1159/000088478

[40] Asadullah, K., Sterry, W. and Volk, H.D. (2003) Interleukin-10 Therapy—Review of a New Approach. Pharmacological Reviews, 55, 241-269.

https://doi.org/10.1124/pr.55.2.4

[41] Wang, L., Liu, J.-Q., Talebian, F., Liu, Z., Yu, L. and Bai, X.-F. (2015) IL-10 Enhances CTL-Mediated Tumor Rejection by Inhibiting Highly Suppressive CD4+T Cells and Promoting CTL Persistence in a Murine Model of Plasmacytoma. Oncoimmunology, 4, e1014232.

[42] Oft, M. (2014) IL-10: Master Switch from Tumor-Promoting Inflammation to Antitumor Immunity. Cancer Immunology Research, 2, 194-199.

https://doi.org/10.1158/2326-6066.CIR-13-0214

[43] Morris, J.C., et al. (2014) Phase I study of GC1008 (Fresolimumab): A Human Anti-Transforming Growth Factor-Beta (TGFβ) Monoclonal Antibody in Patients with Advanced Malignant Melanoma or Renal Cell Carcinoma. PLoS One, 9, e90353.

https://doi.org/10.1371/journal.pone.0090353

[44] Ishida, H., Muchamuel, T., Sakaguchi, S., Andrade, S., Menon, S. and Howard, M. (1994) Continuous Administration of Anti-Interleukin 10 Antibodies Delays Onset of Autoimmunity in NZB/W F1 Mice. The Journal of Experimental Medicine, 179, 305-310.

https://doi.org/10.1084/jem.179.1.305

[45] DePillis, L., Caldwell, T., Sarapata, E. and Williams, H. (2013) Mathematical Modeling of Regulatory T Cell Effects on Renal Cell Carcinoma Treatment. Discret. Contin. Dyn. Syst. - Ser. B, 18, 915-943.

[46] Sakamoto, N., et al. (2015) Phase I Clinical Trial of Autologous NK Cell Therapy Using Novel Expansion Method in Patients with Advanced Digestive Cancer. Journal of Translational Medicine, 13, 277.

https://doi.org/10.1186/s12967-015-0632-8

[47] Kronik, N., Kogan, Y., Vainstein, V. and Agur, Z. (2008) Improving Alloreactive CTL Immunotherapy for Malignant Gliomas Using a Simulation Model of Their Interactive Dynamics. Cancer Immunology, Immunotherapy, 57, 425-439.

https://doi.org/10.1007/s00262-007-0387-z

[48] Rosenberg, S.A. (2014) IL-2: The First Effective Immunotherapy for Human Cancer. The Journal of Immunology, 192, 5451-5458.

https://doi.org/10.4049/jimmunol.1490019

[49] de Pillis, L., et al. (2009) Mathematical Model Creation for Cancer Chemo-Immunotherapy. Computational and Mathematical Methods in Medicine, 10, 165-184.

https://doi.org/10.1080/17486700802216301

[50] Arciero, J.C., Jackson, T.L. and Kirschner, D.E. (2004) A Mathematical Model of Tumor-Immune Evasion and siRNA Treatment. Discret. Contin. Dyn. Syst. - Ser. B, 4, 39-58.

[51] Robertson-Tessi, M., El-Kareh, A. and Goriely, A. (2012) A Mathematical Model of Tumor-Immune Interactions. Journal of Theoretical Biology, 294, 56-73.

https://doi.org/10.1016/j.jtbi.2011.10.027

[52] Bian, X., et al. (2006) Increased Angiogenic Capabilities of Endothelial Cells from Microvessels of Malignant Human Gliomas. International Immunopharmacology, 6, 90-99.

https://doi.org/10.1016/j.intimp.2005.08.004

[53] Hahnfeldt, P., Panigrahy, D., Folkman, J. and Hlatky, L. (1999) Tumor Development under Angiogenic Signaling: A Dynamical Theory of Tumor Growth, Treatment Response, and Postvascular Dormancy. Cancer Research, 59, 4770-4775.

[54] Mabarrack, N.H.E., Turner, N.L. and Mayrhofer, G. (2008) Recent Thymic Origin, Differentiation, and Turnover of Regulatory T Cells. Journal of Leukocyte Biology, 84, 1287-1297.

https://doi.org/10.1189/jlb.0308201

[55] Sun, X., et al. (2016) Mathematical Modeling of Therapy-Induced Cancer Drug Resistance: Connecting Cancer Mechanisms to Population Survival Rates. Scientific Reports, 6, 22498.

https://doi.org/10.1038/srep22498

[56] Cameron, M.A., Davis, A.L., Kostelich, E. and Eikenberry, S. (2009) A Mathematical Model of Angiogenesis in Glioblastoma Multiforme. Arizona State University, Tempe.