Model construction

As explained in detail in the methods section, we have constructed a multi-scale agent-based model of the interactions between macrophages, T cells, and tumor cells, along with the cytokines IL-4 and IFN-γ. The model represents a 2D tissue slice, and cells are constrained to a lattice, with one cell per lattice site. Briefly, macrophages are initially present within the tissue and more are recruited to the tumor site as simulation progresses. They are able to differentiate based on cytokine concentration into either an M1-like or M2-like state. Differentiation is based on an intracellular signaling model, which was then reduced to a simple input/output model using a neural network to relate cytokine concentrations to differentiation. T cells are recruited based on tumor cell death, and act to kill tumor cells. They become fully active upon antigen contact, allowing them to kill tumor cells and secrete IFN-γ. Activation is promoted by M1 macrophages and inhibited by M2 macrophages. Tumor cells proliferate and secrete IL-4. With this model, we predict the effects of three macrophage-based immunotherapies (macrophage depletion, recruitment inhibition, and macrophage reeducation via PI3K inhibition) and how they impact the ability of T cells to remove the tumor. A model schematic along with an example spatial distribution of cells is shown in Fig 1. From the spatial distribution (Fig 1B), which shows a sample simulation without treatment, we see that macrophages surround the tumor, primarily in the M2 state, preventing some T cells from reaching the tumor and inhibiting T cells that are recruited to the tumor site. There are active T cells immediately adjacent to tumor cells; however, there are not enough to eliminate the tumor.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Model schematic and representative simulation result.

(A) Model schematic. T cells secrete IFN-γ, which promotes M1 differentiation. Tumor cells and M2 macrophages secrete IL-4, which promotes M2 differentiation. T cells kill tumor cells. M1 macrophages promote T cell function, while M2 macrophages inhibit it. (B) Representative spatial distribution of cells once the tumor reaches the equilibrium state. Tumor cells (black), T cells (yellow), active T cells (red), M0 macrophages (blue), M2 macrophages (purple). There are no M1 macrophages present at the end of this simulation.

https://doi.org/10.1371/journal.pcbi.1008519.g001

Model behavior without treatment

Prior to simulating macrophage-based interventions, we ran simulations without any treatment to have a reference for comparing the efficacy of different treatment strategies. To understand how various parameters relating to the immune response affect model behavior, we sampled over a wide range of values above and below our base parameter values using Latin Hypercube Sampling (LHS), generating a total of 500 parameter sets. Due to the stochastic behavior of the model, each parameter set was simulated 100 times in order to obtain average behavior. We note that LHS generates a series of parameter sets where each parameter value only appears once. Therefore, the sharp spikes seen in the results are due to each parameter value being simulated once, and not the average response at the parameter value. Despite this, we still see trends for certain parameters. We found that macrophage recruitment rate had the only noticeable effect on the fraction of tumors removed by the immune system. Specifically, higher recruitment rates, and thus more macrophages in the system, lead to a lower fraction of tumors that were removed (S1 Fig). We also found that macrophage recruitment rate and lifespan were correlated to the final number of tumor cells at the end of simulation (S2 Fig). That is, having more macrophages, which would be in the M2 state (discussed in detail below), present in the environment, due to either increased recruitment or longer lifespan, leads to more tumor cells at the end of simulation. These results support the finding that high numbers of TAMs have a worse clinical outcome [15]. This analysis also supports our focus on predicting the effects of macrophage-based therapies. The final values of the model parameters were set so that few tumors were removed by the immune system without treatment, so that tumor killing in subsequent simulations would be due primarily to treatment, and to match those used in similar models.

We used the base model to explore the dynamics of each cell type and diffusible factor present in the TME. Shown in Fig 2 are the time courses for 100 simulations over 200 days without treatment. For all simulations where the immune system failed to eliminate the tumor (99%), each cell type reached an equilibrium state. This state can be considered the “immune control” phase of tumor growth, which is thought to take place over several years and involves the selection of tumor cells that are resistant to the immune system [3]. We conclude that this equilibrium state is indeed due to the cytotoxic function of the T cells and not spatial inhibition, with S3 Fig comparing tumor growth curves in the absence of immune cells and with immune cells present but without function, finding these two curves to be almost identical. After an initial drop due to the introduction of T cells to the environment, the tumor cell population evens out at around 350 cells (Fig 2A). The naive macrophage population drops immediately and stays at a low number due to continuous recruitment to the TME and subsequent differentiation (Fig 2B). In this set of simulations, we see no macrophage differentiation to the M1 phenotype (Fig 2C) whereas the number of M2 macrophages reaches a high level (Fig 2D). The time courses for the total number of T cells (Fig 2E) resembles a delayed version of those for tumor cells. The time courses for active T cells follows this, though at very low levels (Fig 2F). Due to the number of tumor cells, average and maximum IL-4 levels are relatively high (Fig 2G and 2H) while IFN-γ levels (Fig 2I and 2J) are low due to low T cell activation. This difference causes the absence of M0 differentiation into the M1 phenotype. We also briefly examined how initial macrophage density impacts tumor dynamics (S4 Fig). We found that while this impacts initial tumor dynamics, long-term dynamics were not impacted by initial macrophage count.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Time courses with baseline parameters and no treatment.

Time courses for tumors that were not removed by the immune system are shown in black; those for tumors that were removed are shown in red. There is only one simulation here that led to tumor removal. (A) Cancer cells, (B) M0 cells, (C) M1 cells, (D) M2 cells, (E) total T cells, (F) active T cells, (G) average IL-4, (H) maximum IL-4, (I) average IFN-γ, (J) Maximum IFN-γ. (G)–(J): units are number of molecules.

https://doi.org/10.1371/journal.pcbi.1008519.g002

Effects of continuous treatment

After evaluating tumor growth in the absence of treatment, we evaluated the efficacy of three macrophage-based therapies: (1) macrophage depletion, (2) inhibition of macrophage recruitment, and (3) re-education of macrophages via PI3K inhibition. To initially explore how each treatment strategy affects the immune response, we simulated each treatment at ten different strengths. Treatment was started 100 days after tumor initiation, which was after the system had reached equilibrium. The treatment continued until either the tumor was removed or until 200 days of simulation time was reached. We repeated each simulation 100 times and averaged the results to obtain the general effect of each treatment. A primary output of these simulations is the fraction of tumors removed—the number of tumors eliminated after treatment was started. We also calculated the time from start of therapy to tumor elimination and the maximum numbers of M1 macrophages, T cells, and active T cells. Because treatment starts once equilibrium is reached, the system is already at the maximum number of tumor cells and M2 macrophages, so we do not record their maximum values during treatment.

First, we investigated the effects of macrophage depletion (Fig 3A). The “Depletion Probability” displayed on the x-axis is the probability that a macrophage has of being removed at each time step, with an equal probability for each macrophage phenotype. Fig 3A-i shows the percent of tumors eliminated after treatment is started. We see that there is a drastic jump from almost no effectiveness at a probability of 0.001 to almost complete tumor removal at a probability of 0.002, which then increases slightly and stays at complete tumor removal as the probability is increased. It is interesting that there is such a sharp increase in effectiveness without a corresponding increase in the maximum number of M1 macrophages or active T cells. We believe that this increase is due to removal of enough macrophages around the tumor to allow T cells near the tumor to infiltrate better and remove the tumor. Following this, the average number of days needed to remove the tumor decreases and then remains static (Fig 3A-ii). We see that the maximum number of M1 macrophages increases with depletion probability to a point, then begins to decrease at the highest depletion probabilities (Fig 3A-iii). We believe that the increase is due to the decrease in IL-4 as tumor cells and M2 macrophages are removed, allowing newly recruited macrophages to differentiate into the M1 state. In S5 Fig we show the time courses at a depletion probability of 0.006, which shows a decrease in IL-4 and and increase in IFN just prior to the increase in M1 macrophages. As depletion probability further increases, these new macrophages are removed fast enough to decrease the maximum number of M1 macrophages. However, due to the low numbers of M1 macrophages and large error bars, we cannot make any definitive inferences from this. The total number of T cells increases very slightly with depletion probability (Fig 3A-iv) while T cell activation sees a more drastic increase (Fig 3A-v) corresponding with the number of M1 macrophages.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Effects of continuous immunotherapy started at 100 days of simulation.

(A) Macrophage depletion, (B) recruitment inhibition, and (C) PI3K inhibition. (i) fraction of tumors removed after starting therapy. (ii) time (days) from starting treatment to tumor removal. It is averaged over the 100 simulations and is equal to zero if no tumors were removed at that treatment level. (iii) maximum number of M1 macrophages. (iv) maximum number of total T cells. (v) maximum number of active T cells. Note the differences in y-axis scales across treatment strategies. Asterisks signify that a result is statistically significant (p<0.01) from the result of the lowest treatment strength. We note that for the time needed to remove the tumor (ii), we plot the time averaged over only simulations where the tumor was removed. Therefore, while some bars may appear much higher than that of the lowest treatment strength, they only represent a small number of simulations out of 100 and thus were not found to be statistically significant.

https://doi.org/10.1371/journal.pcbi.1008519.g003

Examining the time course for a depletion probability of 0.002 (S6 Fig), we find that in simulations that fail to remove the tumor, the tumor cell population reaches a new equilibrium state that is lower than the original (S6A Fig). As the tumor is removed, there is an increase in the number of naïve macrophages (S6B Fig) due to less of an influence from the tumor for macrophage differentiation. At this level of depletion, no M1 macrophages arise during simulation (S6C Fig), while the number of M2 macrophages decreases (S6D Fig). We also see little difference in T cell and active T cell dynamics (S6E and S6F Fig) between simulations that removed the tumor and those that did not, presumably due to the lack of M1 macrophages, which promote T cell activation.

The second treatment strategy we employed was inhibition of macrophage recruitment, which is similar to macrophage depletion in that it eliminates macrophages regardless of phenotype. On the x-axes for Fig 3B is the “Inhibition Strength,” which is the fraction that the parameter for recruitment rate was reduced by. We see that below 0.6 inhibition, treatment is ineffective in removing the tumor. At 0.6 there is very minor effectiveness, which then rapidly increases with inhibition strength (Fig 3B-i). As inhibition strength increases, we see a slight decrease in the average time needed to remove the tumor (Fig 3B-ii), however it does take longer to remove the tumor than with macrophage depletion probabilities that achieved similar tumor removal (compare Fig 3A-ii and Fig 3B-ii). Whereas macrophage depletion showed some increase in the number of M1 macrophages, they are essentially nonexistent with recruitment inhibition (Fig 3B-iii). We believe this is due to the fact that macrophage depletion removes macrophages from the environment while recruitment inhibition prevents new macrophages from entering. While both methods lower the total number of macrophages, the macrohpages that remain even with recruitment inhibition are closer to the tumor, where IL-4 levels are highest, which keeps them in the M2 state. With macrophage depletion, new macrophages enter the environment farther from the tumor where IL-4 levels are lower, allowing some of them to differentiate to the M1 state. With this, we see that there is no change in the numbers of total or active T cells (Fig 3B-iv,v), with T cell activation being fairly low.

Examining the time courses for a recruitment inhibition of 0.7, which eliminated roughly half of the tumors, we find that the time courses for tumors that were removed are very similar to those that were not (S7 Fig). The dynamics for the tumor cells are very similar for the initial decrease, with tumors that were eliminated continuing to decrease, while tumors that were not suddenly leveling off to a new equilibrium (S7A Fig). We believe this is due to inherent randomness and M2 macrophages preventing T cells from advancing t the tumor. We do see some increase in naive macrophages in tumors that were eliminated (S7B Fig). Like the simulations without treatment, no M1 macrophages appear (S7C Fig). The numbers of M2 macrophages, total T cells, and active T cells exhibit similar behavior as the tumor cells (S7D–S7F Fig). IL-4 and IFN-γ concentrations mirror tumor cell and active T cell dynamics, respectively (S7G–S7J Fig).

The final treatment strategy, PI3K inhibition, yielded similar efficacies as recruitment inhibition, though with a steeper shift from non-effectiveness to effectiveness (Fig 3Ci). However, PI3K inhibition displayed the fastest tumor removal time (Fig 3Cii) due to increases in M1 macropahges and thus T cell function. As intended, because PI3K inhibition aims to reeducate the macrophages, the number of M1 macrophages is much higher with this strategy (Fig 3Ciii). This increase correlates with an increase in the number of total and active T cells (Fig 3Civ,v). We also examined the time courses for a PI3K inhibition of 0.8 (S8 Fig). There is the expected decrease in tumor cells (S8A Fig) and, since this treatment does not remove macrophages, there is an increase in naive macrophages (S8B Fig) due to a decreased pressure to differentiate as the number of cancer cells is reduced. The number of M1 macrophages rapidly increases immediately after treatment, while the number of M2 macrophages rapidly decreases (S8C and S8D Fig). Both total and active T cells increase greatly with the increase in M1 macrophages (S8E and S8F Fig) due to the influence of the M1 macrophages. As with other treatments, IL-4 and IFN-γ levels follow tumor and T cell dynamics (S8G–S8J Fig).

Overall, we find that, at the baseline parameters, treatment efficacy for all three strategies steeply increases from ineffective to very effective when simulated continuously. PI3K inhibition, which converts the macrophages to the M1 state, leads to higher levels of T cell activation and removes the tumor at a faster rate than the other two treatments. Interestingly, although both of these treatments reduce the number of macrophages in the TME, we see that macrophage depletion leads to a slight increase in M1 macrophages while recruitment inhibition does not.

Treatment cycles

While it is ideal to be able to give a continuous treatment, this is unrealistic due to various reasons, such as patient compliance, inability to give treatment continuously due to pharmaceutical implementation, or the desire to give as little treatment as possible in order to avoid potential side-effects. Therefore, we ran several sets of simulations where treatment was cycled on and off. We vary both the total length of the treatment cycle and the amount of time during the cycle that treatment is given. We performed these simulations only for higher treatment strengths, as it can be reasoned that cycling treatment would not be as effective as constant treatment.

For macrophage depletion, we implement treatment by removing the fraction of macrophages equal to “Depletion Strength” at the beginning of each treatment cycle (Fig 4A), shown on the y-axis. The duration of each cycle is shown on the x-axis. Even with cycling treatment on and off, we find that there is a steep increase in the treatment being ineffective to very effective as depletion strength is increased or cycle duration is decreased (Fig 4A-i). The time needed for the treatment to remove the tumor also decreases slightly as depletion strength is increased and cycle duration decreased (Fig 4A-ii). Surprisingly, we see a slight increase in M1 macrophages, total T cells, and active T cells (Fig 4A-iii to 4A-v), especially at the highest depletion strength and slightly longer cycle durations, to a higher level than with constant treatment. This is most likely due to the decrease in IL-4 as tumor cells and M2 macrophages are removed, allowing some of the new macrophages to differentiate into the M1 state.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Effects of cycled immunotherapy started at 100 days of simulation.

For macrophage depletion (A), the fraction of macrophages removed at the beginning of each cycle is given as “Depletion Strength” and the length of each cycle is “Cycle Duration.” For recruitment inhibition (B) and PI3K inhibition (C), the number of days in the cycle that treatment is on for is given as “Days Treatment is On.” Recruitment inhibition is simulated at a strength of 1.0 (complete inhibition) and PI3K inhibition is simulated at a strength of 0.8. For recruitment inhibition and PI3K inhibition, spaces marked with an X are those where treatment-on time is equal or greater to the cycle duration, thus were not simulated. (i) fraction of tumors removed after starting therapy. (ii) time (days) from starting treatment to tumor removal. It is averaged over the 100 simulations and is equal to zero if no tumors were removed at that treatment level. (iii) maximum number of M1 macrophages. (iv) maximum number of total T cells. (v) maximum number of active T cells.

https://doi.org/10.1371/journal.pcbi.1008519.g004

Whereas macrophage depletion continues to be effective when treatment is cycled, recruitment inhibition becomes rather ineffective, even at complete inhibition of macrophage recruitment (Fig 4B). When cycling recruitment inhibition, treatment was turned on for the number of days shown on the y-axis and turned off for the remainder of the cycle (x-axis). When treatment is given for almost the entire cycle, there is strong tumor removal, however this quickly falls when treatment is not given for as long or when cycle duration increases (Fig 4B-i). The time needed to remove the tumor is fairly consistent in cases where there is tumor removal (Fig 4B-ii). While there are little to no M1 macrophages for almost all of the treatment combinations, there is one instance where is a large number of M1 macrophages (Fig 4B-iii). However, because that combination saw very little tumor removal, we assume it is a very rare stochastic occurrence. The numbers of T cells and active T cells remain low and constant across combinations, except the instance where there was an increase in M1 macrophages (Fig 4B-iv,v).

Interestingly, a PI3K inhibition of 0.8, which is the lowest successful inhibition for continuous treatment, is still very successful when cycling treatment, even when treatment is given briefly for long cycle durations (Fig 4C-i). For the bulk of the treatment combinations, the time needed to remove the tumor is fairly constant (Fig 4C-ii). M1 macrophages, T cells, and active T cells are all at their maximum values (Fig 4C-iii to 4C-v). To determine why PI3K inhibition continues to be successful even when given for short durations, we examined the time courses when treatment was given for 2 days out of a 25-day cycle (Fig 5). Interestingly, treatment cycles do not always remove the tumor; however, when tumor elimination does happen, it occurs very swiftly (Fig 5A). It is unclear why this occurs, though we do see that on cycles that do not lead to removal, there is little change in the number of M2 macrophages. Because the sustained response that leads to tumor removal is due to T cell produced IFN-γ, we believe that not enough T cells are becoming activated during these cycles to cause the sustained response. We have highlighted a single time course in red to make it easier to visualize. Immediately following a successful cycle of treatment, we see a gradual rise in naive macrophages as the tumor is eliminated (Fig 5B), and a rapid increase in M1 macrophages followed by a gradual decrease (Fig 5C). M2 macrophages, on the other hand, rapidly decrease since they are being converted in the M1 state (Fig 5D). Corresponding with M1 behavior, the numbers of T cells and active T cells rapidly increase due to influence from the M1 macrophages (Fig 5E and 5F). IL-4 and IFN-γ levels correlate with tumor and T cell dynamics (Fig 5G–5J). We see that due to sustained T cell activation there is sustained IFN-γ, which is responsible for promoting M1 differentiation. This creates a feedback loop between the M1 macrophages and T cells that sustains M1 differentiation and tumor removal even after treatment is removed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. PI3K inhibition of 0.8 for cycled immunotherapy.

Individual time courses for PI3K inhibition of 0.8 at a cycle duration of 25 days with treatment given for 2 days per cycle. Tumors that survived to the end of simulation are shown in black. Tumors that were eliminated are shown in green. One time course is shown in red for ease of understanding. (A) Cancer cells, (B) M0 macrophages, (C) M1 macrophages, (D) M2 macrophages, (E) T cells, (F) Active T cells, (G) Average IL-4, (H) Maximum IL-4, (I) Average IFN-γ, (J) Maximum IFN-γ.

https://doi.org/10.1371/journal.pcbi.1008519.g005

Overall, when cycling treatment on and off, we find that macrophage depletion yields expected results, with higher depletion strengths and shorter treatment cycles leading to a stronger tumor removal. Interestingly, recruitment inhibition becomes very ineffective when cycled, unless treatment is given for almost the entirety of simulation. PI3K inhibition, however, is very effective for almost every combination of time on and time off. We find that this is because converting the macrophages to the M1 phenotype promotes T cell activation and IFN-γ secretion, which sustains the M1 phenotype after treatment is removed.

Changing tumor proliferation rate and macrophage recruitment rate

Of interest is how tumor proliferation rate and macrophage recruitment rate affect the immune response with and without treatment. Presumably these two parameters would lead to a different equilibrium state, which would change the effectiveness of each treatment. We repeated the above analysis for three sets of tumors: increased proliferation rate, increased macrophage recruitment rate, and both. We first consider tumor growth without treatment, and then implement the three macrophage-based treatment strategies.

Compared to the baseline parameters, we see only a slight increase in the equilibrium number of tumor cells (Fig 6A) when tumor proliferation rate is increased from 0.8 per day to 1.2 per day. What is most interesting is that at random points throughout the equilibrium state, the immune system will suddenly remove the tumor, a phenomenon not seen with the baseline parameters. As with previous simulations, a decrease in tumor cell population is followed by an increase in naive macrophage population (Fig 6B). For simulations that removed the tumor, we see a rapid increase in M1 macrophages with a corresponding decrease in M2 macrophages (Fig 6C and 6D). The dynamics of total T cells and active T cells follows that of M1 macrophages (Fig 6E and 6F). IL-4 and IFN-γ dynamics follow tumor cells and T cells, respectively (Fig 6G–6J).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Time courses with increased tumor proliferation rate and no treatment.

The tumor proliferation rate is increased to 1.2/day. Time courses for tumors that were not removed by the immune system are shown in black; those for tumors that were removed are shown in red. (A) Cancer cells, (B) M0 cells, (C) M1 cells, (D) M2 cells, (E) Total T cells, (F) Active T cells, (G) Average IL-4, (H) Maximum IL-4, (I) Average IFN-γ, (J) Maximum IFN-γ.

https://doi.org/10.1371/journal.pcbi.1008519.g006

At the increased tumor proliferation rate, we see a more gradual increase in the effectiveness of constant PI3K inhibition (Fig 7A). However, these results are conflated due to propensity of the immune system to spontaneously remove the tumor during the equilibrium state, as evident by the large error bars in Fig 7B. Across inhibition strengths, there is a large amount of M1 differentiation, T cell numbers, and T cell activation (Fig 7C–7E). Cycling PI3K inhibition is successful throughout the different cycle durations, even with the increased tumor proliferation rate (S9C Fig). Macrophage depletion and recruitment inhibition do not show a clear trend and the results are again conflated due to the inherent tumor removal found without treatment (S9A and S9B Fig, S10 Fig and S11 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Constant PI3K inhibition at tumor proliferation rate of 1.2/day.

(A) fraction of tumors removed after starting therapy, (B) average time needed to remove the tumor, (C) the maximum number of M1 macrophages, (D) the maximum number of total T cells, (E) the maximum number of active T cells. Note the differences in y-axis scales across treatment strategies. Asterisks signify that a result is statistically significant (p<0.01) from the result of the lowest treatment strength. We note that for the time needed to remove the tumor (B), plotted is the time averaged over only simulations where the tumor was removed. Therefore, while some bars may appear much higher than that of the lowest treatment strength, they only represent a small number of simulations out of 100 and thus were not found to be statistically significant.

https://doi.org/10.1371/journal.pcbi.1008519.g007

Whereas increased tumor proliferation only caused a slight increase in the equilibrium tumor population, doubling the rate of macrophage recruitment greatly increased the equilibrium tumor population to roughly twice of the baseline (Fig 8A). Naive and M1 macrophage dynamics (Fig 8B and 8C) are very similar to the baseline while dynamics for the remaining cell populations and cytokines reach a higher equilibrium state (Fig 8D–8J).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Time courses with macrophage recruitment rate doubled and no treatment.

(A) Cancer cells, (B) M0 cells, (C) M1 cells, (D) M2 cells, (E) Total T cells, (F) Active T cells, (G) Average IL-4, (H) Maximum IL-4, (I) Average IFN-γ, (J) Maximum IFN-γ.

https://doi.org/10.1371/journal.pcbi.1008519.g008

At the higher macrophage recruitment rate, macrophage depletion therapy bears similar effectiveness as the baseline recruitment rate (Fig 9A), however effectiveness appears at a depletion probability of 0.003 rather than 0.002. The time needed to remove the tumor is constant across depletion probabilities (Fig 9B) and is a bit higher than at the baseline recruitment rate. There is also a large increase in the number of M1 macrophages (Fig 9C), and the peak at a moderate depletion probability becomes more clear. This response is mirrored with the numbers of T cells and active T cells (Fig 9D and 9E), though with a more gradual tail. The other therapies resemble the baseline case (S12, S13 and S14 Figs).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Constant macrophage depletion started at 100 days with macrophage recruitment rate doubled.

(A) fraction of tumors removed after starting simulation, (B) average time needed to remove the tumor, (C) the maximum number of M1 macrophages, (D) the maximum number of total T cells, (E) the maximum number of active T cells. Note the differences in y-axis scales across treatment strategies. Asterisks signify that a result is statistically significant (p<0.01) from the result of the lowest treatment strength. We note that for the time needed to remove the tumor (B), plotted is the time averaged over only simulations where the tumor was removed. Therefore, while some bars may appear much higher than that of the lowest treatment strength, they only represent a small number of simulations out of 100 and thus were not found to be statistically significant.

https://doi.org/10.1371/journal.pcbi.1008519.g009

Increasing both tumor proliferation rate and macrophage recruitment rate leads to the highest equilibrium tumor population (Fig 10A). The dynamics of the other cells and cytokines behave similarly to the previous simulations without treatment (Fig 10B–10J).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Time courses with increased tumor proliferation and macrophage recruitment rates, without treatment.

The macrophage recruitment rate is doubled, and tumor proliferation rate increased to 1.2/day. Time courses for tumors that were not removed by the immune system are shown in black; those for tumors that were removed are shown in red. (A) Cancer cells, (B) M0 cells, (C) M1 cells, (D) M2 cells, (E) Total T cells, (F) Active T cells, (G) Average IL-4, (H) Maximum IL-4, (I) Average IFN-γ, (J) Maximum IFN-γ.

https://doi.org/10.1371/journal.pcbi.1008519.g010

Interestingly, lower treatment strengths of continuous PI3K inhibition are more effective here than at the baseline parameters, despite the higher equilibrium state of the tumor (Fig 11A), and efficacy increases gradually with inhibition strength. However, at the higher inhibition strengths, treatment is not as effective as it was at the baseline parameters (Fig 3Ci). The time needed to remove the tumor is fairly constant across inhibition strengths (Fig 11B). The numbers of M1 macrophages, total T cells, and active T cells all increase gradually following treatment efficacy, however, due to the wide standard deviations, it does not appear to be too significant of an increase (Fig 11C–11E). Neither continuous macrophage depletion nor recruitment inhibition have a significant ability to remove the tumor (S15 and S16 Figs).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Constant PI3K inhibition at increased tumor proliferation and macrophage recruitment rates.

(A) fraction of tumors removed after starting therapy, (B) average time needed to remove the tumor, (C) the maximum number of M1 macrophages, (D) the maximum number of total T cells, (E) the maximum number of active T cells. Note the differences in y-axis scales across treatment strategies. Asterisks signify that a result is statistically significant (p<0.01) from the result of the lowest treatment strength. We note that for the time needed to remove the tumor (B), plotted is the time averaged over only simulations where the tumor was removed. Therefore, while some bars may appear much higher than that of the lowest treatment strength, they only represent a small number of simulations out of 100 and thus were not found to be statistically significant.

https://doi.org/10.1371/journal.pcbi.1008519.g011

While cycling PI3K inhibition remains similar to the baseline case (S17 Fig), we see a very different behavior with macrophage depletion and recruitment inhibition when tumor proliferation and macrophage recruitment rates are increased, compared to baseline case. In Fig 12 we show the effects of cycling macrophage depletion while the effects of recruitment inhibition, which are very similar, are shown in S18 Fig. Interestingly, moderate depletion strengths and treatment cycle lengths or high depletion strengths at long treatment cycles were the most effective at removing the tumor (Fig 12A). This result is very different from previous simulations, where efficacy correlated to the amount of treatment given. This shows that more moderate treatment could potentially be more effective. While the efficacy here still is not very high, it is higher than continuous treatment for the same tumor parameters. Time needed to remove the tumor is fairly constant for cases where there was tumor removal (Fig 12B). The levels of M1 macrophages and T cells mirror the behavior of tumor elimination (Fig 12C–12E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Cycling macrophage depletion with increased tumor proliferation and macrophage recruitment rates.

The macrophage recruitment rate is doubled, and tumor proliferation rate increased to 1.2/day. (A) fraction of tumors removed after starting simulation, (B) average time needed to remove the tumor, (C) the maximum number of M1 macrophages, (D) the maximum number of total T cells, (E) the maximum number of active T cells.

https://doi.org/10.1371/journal.pcbi.1008519.g012

Overall, we find that the rate of macrophage recruitment is more impactful on the equilibrium state than the tumor proliferation rate. We also note that, surprisingly, increasing the tumor proliferation rate leads to spontaneous tumor removal during the equilibrium state. As can be expected, treatment efficacy decreases as the equilibrium state increases, though PI3K inhibition retained strong efficacy at higher treatment strengths. The most interesting result is from the simulations at increased tumor proliferation and macrophage recruitment. Here, when cycling macrophage depletion and recruitment inhibition, we find that moderate treatments are actually more effective, causing a slight increase in M1 macrophages and tumor removal, however their effectiveness is well below that of PI3K inhibition.

Overview of model

The model represents early tumor growth or a small initial metastasis, and we simulate the tumor in 2D, representing a tissue slice. The model consists of three cell types: T cells, macrophages, and cancer cells, represented as discrete agents. Interactions between cell types and cytokines are shown in Fig 1A. The lattice is a 100x100 grid representing a 1.5x1.5mm tissue slice, with each site being a 15-micrometer square, the size of one cell diameter [40]. As such, only one cell can occupy a site at a time. While different cell types do have different sizes, the “one cell per site” assumption is a necessary limitation of on-lattice models, with many similar models making the same assumption [27,29,30]. Off-lattice models are able to easily account for different cell sizes, however these models are more computationally expensive. Fig 1B displays an example simulation showing the spatial distribution of different cell types. To model the production and diffusion of diffusible factors, we have a layer of partial differential equations (PDEs). Parameters are either taken from previous modeling efforts or set based on experimental observations. Our model does not represent a specific tumor type. Instead, it is meant to examine generalized tumor behavior. We describe the model in detail below and the parameter values are listed in S1 Table, along with supporting references.

Diffusible factors

There are three diffusible factors present: a tumor-secreted factor that activates macrophages and two cytokines, IL-4 and IFN-γ. The first is referred to as the macrophage activation factor, which alerts the macrophages of the tumor and causes them to differentiate. This factor is secreted by the tumor cells and primarily represents high mobility group box 1 protein (HMG-B1) [29]. Also secreted by the tumor, and by M2 macrophages, is IL-4, which is an immunosuppressive cytokine that promotes macrophage differentiation into the M2 tumor-promoting phenotype [41]. The final diffusible factor is IFN-γ which is secreted by activated T cells and is a part of the Th1 response [42]. This immune response promotes macrophage differentiation into the M1 immune-promoting phenotype. IL-4 and IFN-γ were chosen because they are typical pro-tumor and anti-tumor cytokines, respectively. Also, they are used as the inputs to the mechanistic model used to determine macrophage differentiation, described in detail in later sections.

The diffusion and secretion of these factors was modeled using PDEs shown in Eq 1.

(1)

Here, C_i is the concentration of diffusible factor i, k_sec,i is the secretion rate of the factor, and Cells(x,y) are the coordinates of cells that secrete the factor. To solve the equations, we use a finite difference method to discretize them.

Macrophages

Macrophages are initially present in the tissue, and more are recruited to the TME due to the secretion of various chemokines, such as CCL2[15,16,43]. To simulate macrophage infiltration of the tumor, we set a rate of macrophage recruitment [29]. Macrophages enter the TME in the naïve state (M0). Differentiation occurs once a sufficient level of activating factor is present. In our model, a macrophage checks the levels of IL-4 and IFN-γ that are present in their local environment, and then differentiates according to the intracellular model described in the following section. Macrophages migrate towards the tumor, mimicking chemotaxis, and they have a finite lifespan.

Macrophage differentiation model

To increase the level of biological detail in the model, we incorporate a mechanistic ordinary differential equation (ODE) model of macrophage intracellular signaling in response to IL-4 and IFN-γ, which was developed by Zhao et al [41]. Although macrophage phenotype is shown to be on a continuous spectrum [31], we use discrete phenotypes for modeling simplicity, as other papers have done [29,30]. Following the paper by Zhao et al., we first take the model outputs, which are the time courses for iNOS, TNF-α, CXCL9, CXCL10, and IL-12, which are characteristic of an M1 phenotype, and IL-10, Arg-1, and VEGF, which are typical of an M2 phenotype, normalized to their starting values. The product of the time courses of the M1 outputs is divided by the product of the time courses of the M2 outputs to obtain a time course of the “M1/M2 Score.” While there are many possible ways that macrophage differentiation can be determined from here, for our purposes, we calculate this ratio over a 24-hour simulation and then take the average value. If the average value is greater than one, we assume that there is a greater M1-promoting signaling, and the macrophage differentiates into the M1 phenotype. If the value is less than one, the macrophage differentiates into the M2 phenotype. Macrophages are also shown to be plastic, changing their phenotype based on changing environmental conditions [44]. As such, we allow macrophages to reevaluate their local microenvironment, process the input signals via the intracellular signaling model, and redifferentiate every 24 hours. We assume that macrophages can redifferentiate indefinitely.

While ODE models are able to provide an increased level of biological detail to the ABM, they greatly increase the computational burden of simulation. Therefore, to improve computational time, we replaced the mechanistic model with a data-driven model that takes the cytokine concentrations as inputs and outputs the macrophage phenotype. We have already shown that data-driven models are able to predict mechanistic model outputs [45]. Here, we use a neural network to predict macrophage phenotype. Briefly, a neural network is a machine-learning method that predicts outputs given a set of inputs, even with a complex non-linear relationship between the two. It is trained using large sets of input-output data. A neural network can contain several layers, each with many neurons. A single neuron takes as its input the linear combination of the values of either the model inputs or the outputs of the previous layer. The neuron then transforms this input with an “activation function,” such as a sigmoid function. The output of this function becomes the neuron’s output and is used an input to the next layer.

To train the neural network, we ran 100,000 Monte Carlo simulations of the mechanistic model, randomly sampling cytokine concentrations over the range that they would be present in the ABM. We then determined the phenotype as described above, creating a dataset of inputs and outputs. The neural network was then trained in Python using TensorFlow [46]. The final neural network consisted of one hidden layer with four neurons using a sigmoid activation function. We trained the network several times to determine its ability to capture the ODE model outputs, each time randomly splitting the dataset into training and testing sets. With each testing set, the neural network achieved a prediction accuracy of >98%. Because of this very high accuracy in predicting the output of the ODE model, we determined that the neural network is an efficient way to simplify the ODE model into a simple input-output model.

T cells

In order to eliminate the tumor, T cells have to be recruited to the TME. Upon cancer cell death, tumor antigen is brought to the lymph nodes, initiating an immune response. To model this process, we implement a T cell recruitment rate for each time step, calculated using the following equations, which were taken from Gong et al [26]. The rate at which T cells are recruited is then multiplied by the duration of each time step to get the number of T cells that are recruited at that iteration, and then that many T cells are randomly placed in available sites on the lattice, since we assume there is sufficient vascularization for immune cell recruitment.

(2)

(3)

The rate, r(t), of T cell recruitment is calculated using the mutational burden, k_a, which is the extent of tumor cell mutation, the basal recruitment rate, r₁, the neoantigen strength, k_i, and the number of cancer cell deaths over an interval, N_c,death. While antigen characteristics can vary widely across tumors, we leave them constant as it is not the focus of this study. The t_delay parameter accounts for T cell priming and trafficking to the tumor site, leading to a delayed response. The t_window parameter sets a time range centered around the current simulation time minus t_delay from which to accumulate the number of dead cancer cells, which influences the rate of T cell recruitment.

T cells can do one of four actions: migrate towards the tumor (representing chemotaxis), become fully active upon contact with antigen on a cancer cell, proliferate if fully active, and kill a nearby cancer cell. Fully active T cells are also able to secrete IFN-γ to influence M1 polarization. Interactions between T cells and macrophages are described in more detail below.

T cell–macrophage interactions

There are numerous ways that macrophages are able to modulate T cell behavior, including cytokine excretion, antigen presentation, and inhibitory ligand expression. To model T cell-macrophage interactions, we condense these behaviors into having macrophages promote or inhibit T cell activation based on their phenotype [9,15,20,47–49]. M1 macrophages are able to activate neighboring T cells, while M2 macrophages are able to prevent neighboring T cells from becoming active. This yields the following equation that determines the probability of a T cell becoming fully active: (4) where antigenPresence equals 1 if there is a tumor cell or an M1 macrophage present to activate the T cell. numM1 and numM2 are the number of M1 and M2 macrophages, respectively, neighboring the T cell. The parameters k and s are scaling parameters. We formulated this equation, and hand-tuned parameters k and s, to give a reasonable range of probabilities based on neighboring macrophages.

Cancer cells

Cancer cells in the model simply proliferate. After a cancer cell’s internal clock has reached the specified proliferation time, the cancer cell will proliferate, unless there is no room for it to proliferate or it has reached the specified lifespan. If there is not sufficient room for a new cancer cell to be placed, it will become quiescent.

T cell killing of cancer cells

When an active T cell is neighboring a cancer cell, the T cell can recognize the antigen on the cancer cell and begin to kill it. Due to the time it takes to kill a cancer cell, both the T cell and the cancer cell are considered to be engaged and do not undergo other processes for the duration of killing [27]. Once the cancer cell dies, it is removed from simulation, and the T cell is free to continue killing, until it reaches the maximum number of cancer cells that it can kill, at which point the T cell becomes exhausted [27].

Initialization of simulation

Simulation starts with 25 cancer cells placed in the center of the lattice, representing the early growth of a tumor or a micrometastasis. A population of tissue-resident macrophages is randomly scattered throughout the remaining lattice sites. We consider the area surrounding the tumor to be vascularized enough to allow macrophages and T cells to be recruited there.

Model simulation steps and implementation

At each time step, we proceed through the following steps. First, we calculate the diffusion and secretion of diffusible factors for the duration of the time step. Then, we recruit more macrophages to the environment and proceed to iterate through each macrophage in the environment in a random order and allow them to carry out their various functions. We then repeat this for T cells. After this, we iterate through the cancer cells and allow them to proceed with their functions. Lastly, we remove all dead cells from the environment. The model was implemented in C++. The full model is available at: https://github.com/FinleyLabUSC/Early-TME-ABM-PLOS-Comp-Bio.

Treatment

As the focus of this study is on interactions between macrophages and T cells, we examine several macrophage-based therapies to see how they impact the ability of the T cells to remove the tumor. For each therapy, we vary the effectiveness, which was implemented as the fraction that the target parameter was reduced by. Treatments were chosen based on targets in the literature, with the two main strategies being to reduce the number of macrophages in the TME and to reeducate M2 macrophages to an M1 phenotype. The three therapies explored here are macrophage depletion, recruitment inhibition, and PI3K inhibition.

Treatment was implemented in two forms. The first was to simulate treatment continuously for the duration of the simulation. The second was to cycle the treatment, with treatment on for several days, then off for several days.

Recruitment inhibition

As described above, macrophages are recruited to the TME by chemokines such as CCL2. As many macrophages in the TME display an immunosuppressive M2 phenotype, it is thought that preventing macrophage recruitment would allow T cells to better be able to remove the tumor [15,16]. We implement targeting macrophage-recruiting chemokines by reducing the rate at which macrophages are recruited to the environment.

Macrophage depletion

Similar to recruitment inhibition, depletion of macrophages can be used to reduce the number of macrophages in the TME [15,16]. It is to be noted that both of these strategies fail to discriminate between M1 and M2 macrophages, thus also inhibiting the immune-promoting properties of M1 macrophages. We model macrophage depletion by giving the macrophages a probability of undergoing apoptosis at each timestep that treatment is on. When simulating continuous treatment, we decreased the range of depletion probabilities until we reached the point where the lowest probability fails to remove a large number of tumors so that we could compare ineffective levels of depletion to effective levels.

Macrophage reeducation

In order to prevent the elimination of immune-promoting M1 macrophages, reeducation of M2 macrophages into an M1 phenotype is a potential therapy. One noted target is inhibition of PI3K[15,16]. As this is present in the mechanistic model that was utilized for macrophage differentiation, we include variations in PI3K activity in the Monte Carlo simulations and have PI3K activity as an input to the neural network. Targeting this parameter allows us to redifferentiate the macrophages into an M1 phenotype. We simulate this treatment by reducing the value of the PI3K activity parameter.