Skip to main content

Modeling cancer immunoediting in tumor microenvironment with system characterization through the ising-model Hamiltonian


Background and objective

Cancer Immunoediting (CI) describes the cellular-level interaction between tumor cells and the Immune System (IS) that takes place in the Tumor Micro-Environment (TME). CI is a highly dynamic and complex process comprising three distinct phases (Elimination, Equilibrium and Escape) wherein the IS can both protect against cancer development as well as, over time, promote the appearance of tumors with reduced immunogenicity. Herein we present an agent-based model for the simulation of CI in the TME, with the objective of promoting the understanding of this process.


Our model includes agents for tumor cells and for elements of the IS. The actions of these agents are governed by probabilistic rules, and agent recruitment (including cancer growth) is modeled via logistic functions. The system is formalized as an analogue of the Ising model from statistical mechanics to facilitate its analysis. The model was implemented in the Netlogo modeling environment and simulations were performed to verify, illustrate and characterize its operation.


A main result from our simulations is the generation of emergent behavior in silico that is very difficult to observe directly in vivo or even in vitro. Our model is capable of generating the three phases of CI; it requires only a couple of control parameters and is robust to these. We demonstrate how our simulated system can be characterized through the Ising-model energy function, or Hamiltonian, which captures the “energy” involved in the interaction between agents and presents it in clear and distinct patterns for the different phases of CI.


The presented model is very flexible and robust, captures well the behaviors of the target system and can be easily extended to incorporate more variables such as those pertaining to different anti-cancer therapies. System characterization via the Ising-model Hamiltonian is a novel and powerful tool for a better understanding of CI and the development of more effective treatments. Since data of CI at the cellular level is very hard to procure, our hope is that tools such as this may be adopted to shed light on CI and related developing theories.

Peer Review reports


The scientific understanding of the interaction between the Immune System (IS) and cancer has come a long way since the formal postulation of the cancer immunosurveillance hypothesis by Burnet & Thomas in the late 1950s. Said development, however, has not been without difficulties due to both technological limitations as well as justified resistance from the scientific community to embrace certain ideas that were only partially supported by the evidence available at the time. By the beginning of the twenty-first century, several important discoveries had guided the evolution of the immunosurveillance hypothesis into the concept of cancer immunoediting (CI); a brief review of that evolution is found in [1, 2]. Two decades later, CI is a well established idea describing the process whereby the IS can both protect against cancer development as well as, over time, promote the appearance of tumors with reduced immunogenicity [3].

Cancer immunoediting comprises three phases that take place in a Tumor Micro-Environment (TME): (1) Elimination (immunosurveillance), wherein the IS detects and largely eliminates abnormal cells; (2) Equilibrium, describing a process (that may occur over several years) during which the IS maintains a selection pressure (in an evolutionary sense) on tumor cells, which in combination with increasing genetic instability (mutations that generate cell heterogeneity), promotes the appearance and survival of tumor cell variants with augmented capabilities to avoid or suppress immune actions (detection and attack); and (3) Escape, wherein the immunologically adapted tumor gradually expands without control until it becomes clinically apparent [4].

Although conceptually simple, the different processes referred to as CI are extremely complex and dynamic, since the interaction between the IS and tumor cells involves numerous components from both the Innate and the Adaptive IS, and the behavior of these components, as well as that of the tumor cells, changes over time. From a practical point of view, the Escape phase, or rather the transition between the Equilibrium phase and the Escape phase, is of most interest because it can be considered as the point at which the IS losses effectiveness and cancer in a strict sense (a disordered replication of malignant cells in an organ tissue [5]), develops. Because of this, current immunotherapy research is focused on acquiring a better understanding of how cancer cells evade or suppress the immune response against them and on finding solutions for this tumor resistance.

Cancer immunoediting

The interaction between Cancer Cells (CCs) and the Immune System (IS) taking place in the primary TME is quite intricate, since it involves numerous kinds of cells (see Table 1), as well as a complex exchange of chemical signals (cytokines, see Table 2) that govern their behavior [6]. Herein we provide a very brief and basic explanation of the CI stages, merely enough to serve as a technical underpinning for our model. Detailed overviews of CI can be found in [1, 2, 7]. The IS is composed of two subsystems [8]: the Innate IS (IIS) and the Adaptive IS (AIS). The IIS comprises the genetically inherited immune actions and includes cell types such as the NKs, M\(\upvarphi\)s and Ns that act as first responders against aberrant cells [9, 10]. The AIS evolves via a training process based on the successful immune responses to previously unknown diseases; it includes T and B cells of many varieties [11]. Other cell types that are involved in CI are DCs and, importantly, MDSCs [12, 13].

Table 1 Immune system elements in the TME
Table 2 Cytokines present in the TME involved in Cancer Immunoediting


Normal cells are transformed into CCs because of carcinogens or mutations and express tumor antigens (NKG2D ligands, MHC-I molecules, etc.) which are recognized by NKs and the cytotoxic CD8+T cells. NKs induce apoptosis in cancer cells via cytotoxic molecules (such as perforin) or antibodies, and recruit M\(\upvarphi\)s and Ns to dispose of dying tumor cells. DCs act as messengers between the IIS and the AIS by processing antigen material and presenting it to the elements of the AIS: CD4+T, CD8+T, NKT, and B cells. Activated T cells and NKs secrete IFN-\(\upgamma\) (inflammatory) which activates M\(\upvarphi\)s and Ns, and inhibits angiogenesis. The activated B cells produce tumor specific antibodies while CD8+T cells induce apoptosis via cytokines and interaction with TRAIL receptors on CCs. \(\upgamma \updelta \,\)T, CD8+T and NK cells use NKG2D receptors to recognize induced-self antigens [6, 14]. In an immunocompetent host with a healthy IS, immunosurveillance continuously eliminates abnormal cells throughout the body and prevents the outgrowth of cancer cells. Nevertheless, some CCs can survive this immune control and enter Equilibrium.


This is an immune-mediated state of quiescence wherein proliferation and elimination rates of CCs equal each other. This state also implies a balance between the production of anti-tumor (IL-12, IFN- \(\upgamma\)) and pro-tumor (TGF-\(\upbeta\) [15, 16], IDO) cytokines. The IS continued response is mostly carried out by CD4+T and CD8+T cells, with participation of NKs and T\(_{\mathrm{reg}}\) cells. Most importantly, during the Equilibrium phase tumor cells undergo an evolutionary process (often described as editing or sculpting) in which the IS activities constitute a form of environmental pressure and mutations act as adaptations to said environment (i.e. Darwinian selection leading to survival of the fittest). The Equilibrium phase may occur throughout several years during which poorly/non-immunogenic and immunosuppresive transformed cell variants emerge [7, 17, 18]. At this point the IS losses effectiveness and the Escape phase begins.


This phase is characterized by tumor growth without hindrance from the IS. For this to occur, a number of complementary and intricate mechanisms are set in motion that essentially causes the TME to function in a diametrically different manner to what it would normally do. CCs implement strategies to: (1) Evade immunorecognition by decreasing their expression of antigens and increasing that of anti-apoptotic molecules; (2) Suppress immune function by secreting cytokines such as TGF-\(\upbeta\), IL-6, IL-10 and NKG2D ligands (that inhibit the cytotoxicity of \(\upgamma \updelta\)T cells); and (3) Recruit immune cells to indirectly promote tumor growth (also by suppressing immune function). Recruited cell types include B, B\(_{\mathrm{reg}}\), NKT, NK, \(\upgamma \updelta\)T and MDSCs (that suppress the activation of T and NK cells) [4, 19].

Scope and contributions

In previous paragraphs, the phases of CI have been conceptualized as dynamic scenarios wherein two rivaling forces face each other. In the Elimination phase the IS prevails over tumor cells; during the Equilibrium phase the opposing forces remain balanced; in the Escape phase tumor cells triumph over the IS. This is of particular significance for the present work because herein we propose to characterize the IS-tumor interaction by means of an energy-based model developed for the study of phase transitions. Specifically, we employ the Ising-model Hamiltonian [20] which we have used before for the analysis of complex interactions between opposing agents [21, 22] and apply it to the analysis of CI. To the best of our knowledge, this is the first model of CI that employs this energy-based approach. Note that the “energy” to which we refer is a unitless measure used to characterize computational modelsFootnote 1 not to be understood as a physical quantity.

Our proposal is developed in the context of systems biology [8, 23, 24] as an Agent-Based Model (ABM) where emergent behavior is driven by Gompertzian growth and the probabilistic interactions between the agents. Just as any model describing the intricate IS-tumor interaction, wherein many processes are not completely understood yet, the present proposal introduces simplifications in order to make the model manageable and understandable: we disregard the tumor physical structure which gives place to further categorization of the CCs according to their location inside a solid tumor; also, this work is limited to the avascular phase of tumor growth. Such aspects have been discussed elsewhere [25,26,27,28,29]. We emphasize that the main objective of modeling and simulation of CI is to promote the scientific understanding of this complex process, by the generation of emergent behavior in silico that is very difficult to observe directly in vivo or even in vitro. In the following section we offer a brief recount of relevant prior efforts towards said objective.

Related work

The tumor-immune system interaction in the TME can be formalized as a biological system [23]. Computational models of such a system are of great importance, since the analysis made possible by simulations can improve our comprehension of the process behind the disease [21, 30, 31].

In the scientific literature on cancer tumor growth [32], as well as on that related to the IS response [33, 34], the system biology modeling is carried out under either the continuous [35] or the ABM approaches [25, 36]. The continuous approaches use ordinary differential equations (ODE) [31, 37] or partial differential equations (PDE) [38, 39] to establish the parameters of tumor growth and control growth dynamics; however, frequently the required initial and border conditions are not well defined, and changes that can emerge during the process are difficult to incorporate [26]; these issues constitute main drawbacks of this approach.

In the ABM approach, simple agents represent the elements of cancer and of the IS [40]; the interaction between these results in a complex emergent behavior which is not restricted to the linear combination of the individual elements [41, 42]; hence ABMs can incorporate the changes that bring about non-deterministic interactive processes [21, 25, 36].

In 2006, Novozhilov et al. [35] described a deterministic model for oncolytic viruses used as anti-cancer therapy. The Lotka-Volterra logistic equations were employed to model the spread of the virus infection in the tumor, and the high proliferation of CCs followed an exponential growth-function at early stages of tumorigenesis. Using the model of the virus therapy it was concluded that both the infected and uninfected tumor cells could be eliminated over time, even to complete recovery.

Studies of in-vitro stimulation of T cells for patient treatment have addressed the loss of immunocompetence with age, particularly to fight cancer. Figueredo and Aickelin [34] employed an ABM simulation to show that the processes of immune system aging causes the populations of naïve T cells (those able to respond to novel pathogens) to decay over time; aging affects the naïve T cells response against cancer as well as the response to the anti-tumour vaccination process.

In 2012, Wilson and Levy [37] reported experiments that suggest that TGF-\(\upbeta\) inhibition could amplify the anti-tumor immune response when combined with a tumor vaccine. The ODE-based model of cooperative interaction follows the dynamics of the tumor size, TGF-\(\upbeta\) concentration, cytotoxic and T\(_{\mathrm{reg}}\) cells.

In 2015 Wells et al. [43] developed a hybrid discrete-continuous computational model of a nascent metastatic tumor to investigate how functional and spatial heterogeneity of cell types impact tumor pathogenesis. They discovered that tumor escape was enhanced by heterogeneity in the responses of individual immune cells to their environment. The authors carried out simulations assuming deterministic or stochastic polarization of M\(\varphi\)s. Stochastic polarization was modeled byFootnote 2: \(p(M_2)=0.5(1+erf((-1+MS2/p13)/p16))\) where erf is the error function, p13 is a threshold on MS2 (an effector cytokine) and p16 is a stochasticity parameter; in turn, \(p(M_1)=1-p(M_2)\). Although the model of Wells et al. is limited to very specific mechanisms (such as macrophage functional polarization occurring within 5 days of tumor implantation), their work is noteworthy because it presented multi-parametric sensitivity analyses through which the capabilities of ABMs to build an understanding of the phenomena being modeled are illustrated.

Ku-Carrillo et al. [44] described a model of cancer tumor growth that includes the IS response, which is weakened by the effects of obesity and highly-caloric diets. Their model employs ODEs and logistic growth of the amount of fat stored in the adipocyte cells of the patients.

In 2017, Pourhasanzade et al. [25] presented an ABM of tumor growth with some resemblance to this work. Their study discussed the immune response to prevent cancer development, focusing on the growth of solid tumors; in contrast, our proposal models more of the IS elements and employs an energy-function to evaluate the state of the system. Thus, these two works could be considered complementary.

In 2019 Norton et al. [6] summarized the applications of ABM and hybrid modeling to the TME and cancer immune response, with an emphasis on intra-tumor heterogeneity and interactions between cancer cells and stromal cells, (including immune cells), on the tumor-associated vasculature in relation to the immune response, and on cancer immunotherapy. In their review, the approaches to model the immune system are broadly categorized into top-down and bottom-up models; top-down approaches include those based on ODE, PDE and stochastic differential equations (SDE), and model population of cells and their mean behavior at the macroscopic level. In contrast, the bottom-up approaches track individual cells (or other microscopic elements) and their interactions, from which complex emergent behaviors arise. Although features such as stochastic behavior and heterogeneity are easier to capture via bottom-up models, these require more computational power to track the individual agents, and thus the number of agents that can be modeled is a function of the computational resources available.

Nuñez-López et al. [45] studied the interplay between tumor cells and the IS starting from a deterministic model and transforming it into an SDE model whose simulations were related to the phases of CI. Previously, other works had compared models based on SDEs, ODEs and PDEs to ABMs, notably Figueredo et al. [46, 47]. In these works, tumor proliferation and death were defined by \(p(T)=aT^\alpha\) and \(d(T)=bT^\beta\), respectively, where T represents the number of tumor cells. The transition rates for proliferation and death are given by the logistic growth \(aT(1-Tb)\) and the authors described the parameters and behaviors of agents involved in tumor growth as follows:


Reactive behavior

Proactive behavior

a, \(\alpha\), b and \(\beta\)

Dies if rate \(< 0\)

Proliferates if rate \(> 0\)

where a, b, \(\alpha\) and \(\beta\) are parameters that can vary substantially, depending on the case study. A tumor cell agent thus possessed two possible states, alive or dead, and would move between states through two possible actions, dies or proliferates (clearly, if an agent moves to state dead it will remain in that state). As can be seen, one advantage of ABMs is the simplicity of the agents and of the imbedded rules for individual behavior.

Very recently, Sajid et al. [48] studied cancer niche construction by means of a 2-D Ising model to capture the interaction between clusters of cancer cells and of healthy cells, employing the Kikuchi free energy approximation to describe the spatiotemporal evolution of said model. Without pointing out a direct connection, their study describes the building blocks required for CI: apoptosis (or cell death, required in Elimination), local cancer growth (a state akin to Equilibrium since cancer is kept local), and metastasis (which occurs in the Escape phase). However, Sajid et al. do not discuss CI itself, their model is not an ABM and they employ a different formulation of the system Hamiltonian, more adequate for their purposes, instead of the mean field approximation [21] used in our study.

The computational model presented herein is a hybrid proposal between ABM and the continuous approach, since the behavior of the different agents in the simulated system is governed by relationships derived in both the continuous and the discrete domains, including a logistic growth function (a main choice in the ODE- or PDE-based modeling) for the cancer and IS elements, and the Ising-model Hamiltonian [20] to characterize the interaction of agents as cooperation among partners or confrontation between opponents [21, 22]. We conclude by pointing out that the works of Torquato [49], Barradas et al. [21], Sajid et al. [48], and the present study can be seen as complementary to each other, and together describe the basis for an energy-based computational model of cancer from the immunological perspective.

Agent-based model of cancer immunoediting

In this work, an Ising model energy function is employed to characterize the interaction between cancer-cells and IS-cells within an agent-based framework. A good model should allow us to generate and observe the different phases of Cancer Immunoediting, while avoiding an inordinate amount of control parameters. Our model requires only the setting of four hyperparameters; after that, the model dynamics and the end result depend entirely on the stochastic behavior occurring in the simulator according to the three main components described below.

Modeling cancer immunoediting

We study a simplified version of the cancer versus IS interaction in the TME, for which the biological system modeled is depicted in Fig. 1 and the corresponding elements are listed in Table 1. In our ABM there are agents representing these elements, with associated rules that govern their behavior. The simulation of the system is carried out by sequentially and iteratively invoking the different types of agents, according to the pseudocode in Algorithm 1. The actions of agents are probabilistically defined, and the complex behavior of the system which we want to characterize arises from their interaction with other agents of different types. Said interaction becomes intricate by the emergence of pro-tumor actions from IS elements. For instance, added complexity occurs when M\(\varphi\)s and Ns are “recruited” to inhibit immune functions (by release of inhibitory cytokines), becoming pro-tumoral [11].

Fig. 1
figure 1

IS-Cancer interaction in TME: some Ns and M\(\upvarphi\)s can become pro-tumoral. T\(_{reg}\) cells regulate CD8+T and CD4+T cells. NKs attack CCs and promote T cells

Although it has been simplified, our model contains numerous parameters, the most important of which are listed in Table 3. Since adjusting each of these parameters individually is undesirable, our approach is to define a small set of four hyperparameters, based on which the values of all the parameters can be set automatically. Internally, all of the parameters are given individual values, but to the user this process is transparent. A detailed explanation of parameter initialization is provided in the “Simulation methodology” section.

Table 3 Model parameters of interaction

The probabilities referred to in this subsection regulate the behavior of the elements involved in the TME. The actions to be carried out by the cancer and IS elements are summarized in three diagrams. A complete system flowchart is shown in Fig. 2; the cancer-IIS interaction is summarized in Fig. 3, and the cancer-AIS interaction is described in Fig. 4.

Fig. 2
figure 2

System flowchart. Points numbered 1 to 4 indicate entry and exit points to the diagrams in Figs. 3 and 4. Blue and red indicate mostly anti-tumoral and pro-tumoral processes, respectively

Fig. 3
figure 3

The cancer growth and the IIS response at the micro-environment

Fig. 4
figure 4

The tumor immune interaction in the TME with the AIS elements CD8+T, CD4+T and T\(_{\mathrm{reg}}\) cells

For every possible action to occur, certain conditions must be met, which are defined in terms of the corresponding probabilities. Over time, and as a result of the different actions, the number of CCs and IS cells will change, and with them, the energy of the system, characterized by the Ising-like Hamiltonian (see “Characterizing the stochastic interaction” section).

If the CCs are not eliminated by the anti-tumor elements of the IS, the rate of tumor growth will depend on the strength of the tumor (expressed through the TGF-\(\upbeta\)) and on the Ns and M\(\upvarphi\)s that become pro-tumoral. Meanwhile, the strength of the IS response depends on the age and health of its elements (NKs, Ns and M\(\upvarphi\)s) and those of the AIS (CD4+T, CD8+T, and T\(_{reg}\) cells). If not killed by other agents, the active period of the cells will be probabilistically adjusted based on their age.

It has been suggested that the probability of successful interaction between CCs and the anti-tumor neutrophils N\(_1\) is complementary to the probability of interaction between tumor cells and the pro-tumor neutrophils N\(_2\) [11]. A similar relationship holds for the probability of interaction between the tumor cells and the anti-tumor macrophages M\(_1\) and the complementary probability of interaction between the tumor cells and the pro-tumor macrophages M\(_2\).

figure a

Modeling growth in cancer immunoediting

A needed complementary component in the above description is a model for the tumor growth and for the recruitment of IS elements, which define the strength of the cancer and of the immune response, respectively. We describe the growth of the cell populations by means of logistic functions, which approximate Gompertzian growth, a popular choice in modeling growth patterns, in particular for modeling the population of cells in tumor growth [39, 50,51,52]. The logistic function f(t) with midpoint \(t_0\), maximum value of a, and growth rate k, is given by:

$$\begin{aligned} f(t) = \frac{a}{1 + e^{-k(t - t_0)}} \end{aligned}$$

which has \(\mathbb {R}\) as the usual domain, but as time is positive, we translate it to \(\mathbb {R}^{0+}\) (see the bottom right plot in Fig. 5).

Fig. 5
figure 5

Interface of our simulator with results of one simulation illustrating the Equilibrium phase

Having defined the growth of a cluster of cells as f(t), the growth rate is given by the derivative of said function; letting \(\varphi (t)=1/(1+e^{-k(t-t_0)})\) we have:

$$\begin{aligned} f'(t) = a\varphi '(t)=\frac{-a\left( 1+e^{-k(t-t_0)} \right) '}{\left( 1+e^{-k(t-t_0)} \right) ^2}= \frac{ak \left( e^{-k(t-t_0)}\right) }{\left( 1+e^{-k(t-t_0)}\right) ^2} = ak\varphi (t)(1-\varphi (t)) \end{aligned}$$

Through \(\gamma (t):= f'(t)\) (a logistic distribution), the size of the cell clusters can be updated over time by means of the expression:

$$\begin{aligned} n(t + \Delta t) = \gamma (t + \Delta t)n(t) \end{aligned}$$

where \(\Delta t\) represents a positive interval of time of arbitrary length, and the domain of n(t) is also \(\mathbb {R}^{0+}\). The cell clusters are revisited in (5) and (6) below.

Modeling actions and interactions of agents

To model the actions of agents listed in Algorithm 1, first let us distinguish between actions involving two agents, called interactions, and actions involving only one agent, like cell aging or one anti-tumoral cell becoming pro-tumoral. For all interactions, an agent is invoked and paired with a probabilistically-chosen second agent of adequate target-type (for instance, a M\(\varphi\) is paired with a CC) and their interaction (with anti- or pro- tumoral effect depending on the polarization of the M\(\varphi\) in this example) occurs according to the following rule:

$$\begin{aligned} \text{Interaction} = {\left\{ \begin{array}{ll} True, &{} p_k \ge r, r \in [0, 1]\,\text{is a uniform random variable.} \\ False, &{} \text{otherwise.} \end{array}\right. } \end{aligned}$$

where \(p_k\) represents the applicable probability from those listed in Table 3. Equation (4) describes a generic rule for interactions. The particular interaction rules are described in Table 4, among the rules for single-agent actions.

Table 4 Rules for actions and interactions between agents (\({\mathcal {N}}\): Normal distribution)

Aging of the IS cells leads to loss of efficiency against cancer [34]. In this work, individual cell age is modeled by dividing a cell’s life into three parts: the first third includes the time from cell division to the beginning of maturity; the second third comprises the period of time during which the cell is fully developed and works most efficiently; the last third includes the cell’s aging and death; T\(_{reg}\) cells have a role in regulating or suppressing other cells in the IS, such as CD4+T and CD8+T cells [53]. Below we describe the way in which these interactions are characterized.

Characterizing the stochastic interaction

The Ising model is a classic in the formalization of ferromagnetism that arises from the interaction between molecule spins [14]; it is a simple model from statistical physics used to describe energy interactions and phase transitions of matter. The 2-D Ising model Hamiltonian H, given in (5), is increasingly employed in the modeling of complex processes in biology, chemistry and their interdisciplinary matters [49, 54]. In our present context, \(x_i\) and \(x_j\) represent two interacting clusters of cells, belonging to the tumor and to the IS. The Hamiltonian quantifies the energy involved in the interaction between those clusters, which determines the outcome of the interaction, i.e. whether the CCs or the IS cells dominate the interaction.

$$\begin{aligned} H = -\frac{1}{2}\sum _{i,j}^{n}w_{ij}x_{i}x_{j} - \nu \sum _{i}^{n}h_{i}x_{i} \end{aligned}$$

The characterization of each cluster of cells is given by:

$$\begin{aligned} x_i = c_i n_i \end{aligned}$$

where \(c_i \in \{1, -1\}\) according to whether the cluster \(x_i\) acts in favor of the IS or the cancer, respectively. The variable \(n_i\), computed by (3) is the number of cells in \(x_i\); notice that this may refer to groups of cells or to a whole population of them.

The energy in (5) is also determined by variables \(w_{ij}\), which weight the interaction of cells that cooperate or antagonize each other as members of \(x_i\) or \(x_j\). Finally, \(\nu\) denotes the whole (magnetic in ferromagnetism) field energy, and \(h_i\) the way in which this affects each cluster of cells \(x_i\). In the present application, these parameters are related to the effect of the combined chemical signals in the micro-environment. The specific values for the weights and the field energy require to be specifically tuned to experimental data, which for simplicity is not used in this work. Instead we set \(w_{ij}=1\) and \(\nu =0\), which will allow us to analyze the unweighted and unbiased interaction of the elements in our system.

The model Hamiltonian (5) is a suitable measure of the model energy, which indicates the dominating side in a system of opposing forces at any time during the simulation.

Simulation methodology

With the purpose of demonstrating the potential uses of our proposed model, we carried out simulations to illustrate how the model can be employed to analyze the interaction between cancer cells and the immune system through time. In particular, we are interested in showing that the model is capable of generating the different phases of Cancer Immunoediting, that the phases occur for reasonable values of the model’s hyperparameters, and that these phases can be correctly characterized by the model’s Hamiltonian. The model was implemented in the NetLogo multi-agent modeling environment [55] and includes a graphical interface through which a user can input the model’s hyperparameters and observe the outcome of one individual simulation as a set of plots (Fig. 5) that show (a) the number of anti-tumoral and pro-tumoral cells through time; (b) the number of each type of cell in the CI system; (c) the Hamiltonians for the anti-tumoral and pro-tumoral portions as well as the total Hamiltonian in the TME; (d) the tumor growth rate.

The four hyperparameters, two for the IS \((\mu _1, \sigma _1)\) and two for the CCs \((\mu _2, \sigma _2)\) can be introduced into the simulator by means of slider controls; these variables can be set to any values between 0 and 1 to define two Normal probability distributions. Next, scaling factors \(\eta _k\in [0, 1]\) are sampled from the distributions defined by \((\mu _1, \sigma _1)\) or by \((\mu _2, \sigma _2)\) depending on whether the k-th parameter (of those listed in Table 3) relates to the IS cells or the tumor cells, respectively. In this way, the individual probabilities and strength levels of the elements of the IS and the cancer can be adjusted automatically, by setting:

$$\begin{aligned} \text{Parameter}_k = \eta _k\cdot b_k \end{aligned}$$

where \(b_k\) is the base value for \(\text{Parameter}_k\).

Notice that this procedure is not part of the model, but a methodology to simulate the model under different sets of parameters to carry out the system characterization. A similar procedure was followed in [43] to characterize the TME network robustness and the role of heterogeneities, but it was done by sampling the values of their parameters through a linear distribution. In this study it is more convenient to sample the parameter values under Normal distributions centered around reference values chosen for each of the phases of CI, because this facilitates the organization of the simulations and allows us to control the variance of initial parameters between simulations; specifically, for all the simulations we fixed \(\sigma _1= \sigma _2=0.05\). Keeping the variance of the initial parameters fixed becomes relevant for the analysis of the results, because then the observed differences in variability between CI phases can be attributed to the functioning of the model, rather than to the way in which the initial parameters were distributed.

As stated, herein we are particularly interested in generating the scenarios corresponding to the three phases of Cancer Immunoediting: Elimination, Equilibrium and Escape. For the Elimination phase the strength of the IS was set significantly higher than that of the cancer (\(\mu _1 = 0.7\), \(\mu _2 = 0.3\)) since the objective is to simulate an IS operating optimally by consistently eliminating a moderate amount of abnormal cells. For the Equilibrium phase the strength of both subsystems was lowered and set much close to each other, with that of the IS still being higher than that of cancer (\(\mu _1 = 0.3\), \(\mu _2 = 0.15\)); this corresponds to an IS that keeps in check, but just barely, a small number of cancer cells. For the Escape phase the strength of the cancer was increased, but most importantly, the strength of the IS was set to a value for which the IS remains operational although it cannot prevent the survival and reproduction of CCs (\(\mu _1 = 0.4\), \(\mu _2 = 0.6\)). The distributions are shown in Fig. 6.

Fig. 6
figure 6

Hyperparameter distributions used to model the three phases of CI

The simulations described above are designed to show the robustness of our model to variability of its parameters around central values (this variability is called Variance). In a second batch of simulations the central values (the means of the distributions used to sample the hyperparameters) were displaced to show that our model is also robust to changes of the central values (this is described as Bias). Overall, our model was tested on 15 sets of hyperparameters, and for each set of hyperparameters 50 simulations were performed to provide statistical support.


The results presented in this section come from a probabilistic system and thus a meaningful analysis must consider the average result and dispersion over many simulations. Representative results (over 50 simulations per scenario) regarding the number of cells and the Hamiltonians are shown in Figs. 7, 8 and 9 for the Elimination phase, Equilibrium phase and Escape phase, respectively. The shaded region around each curve illustrates the standard deviation of the results. The sign of the Hamiltonian is positive for the pro-tumoral agents in the system and negative for the anti-tumoral agents (signs are arbitrarily assigned and can be reversed with no effect on the final conclusion). The figures also show the sum of the two Hamiltonians: values close to zero indicate that the forces of the adversaries are nearly balanced at that particular time.

Fig. 7
figure 7

Simulation of the elimination phase of cancer immunoediting

Fig. 8
figure 8

Simulation of the equilibrium phase of cancer immunoediting

Fig. 9
figure 9

Simulation of the escape phase of cancer immunoediting

Figures 10, 11 and 12 illustrate the effect of using different hyperparameters on the number of cells and the Hamiltonians while simulating the Elimination, Equilibrium, and Escape phases, respectively. The objective is to observe how sensitive is our model to the setting of its control parameters.

Fig. 10
figure 10

Effect of hyperparameters variations on simulation results for the elimination phase (\(\mu _1 = \mu _2 + 0.40\))

Fig. 11
figure 11

Effect of hyperparameters variations on simulation results for the equilibrium phase (\(\mu _1 = \mu _2 + 0.15\))

Fig. 12
figure 12

Effect of hyperparameters variations on simulation results for the escape phase (\(\mu _1 = \mu _2 - 0.20\))

Finally, Fig. 13 shows the probability density estimates of the energies in the system corresponding to the second half of the simulations (once the system has converged), presented as violin plots for each of the three CI phases. These density estimates were obtained via kernel density estimation [56]. Table 5 reports the means and medians of the probability densities.

Fig. 13
figure 13

Energy probability density estimates in different phase of CI

Table 5 Central tendency statistics of probability densities: the difference between the means and the medians make the distributions in Fig. 13 left- or right-skewed


Based on the results shown in Figs. 7, 8 and 9 the following observations can be made. First, with respect to the number of anti-tumoral and pro-tumoral cells, it should be noticed that the simulations start with an initial population of these cells and it is only during the first fourth of the simulations that the growth rate of CCs increases; after that, following Gompertzian growth, it decreases gradually. Thus, as time passes, the number of cells naturally tends to diminish due to apoptosis and, in the case of CCs due to elimination by cytotoxic cells (except in Escape).

During the Elimination phase (Fig. 7), an anti-tumoral strength of 0.7 versus a pro-tumoral strength of 0.3 translates into initial cell populations of about 180 anti-tumor cells versus 60 pro-tumor cells. Cell population growth is not observable, because the significantly larger number of anti-tumor cells prevents the proliferation of CCs and in turn the decreasing number of these impedes recruitment of more anti-tumor cells. Notice how the corresponding Hamiltonians reflect this behavior: the energy of the pro-tumoral agents is initially about one third of the anti-tumoral energy, then the pro-tumoral energy falls rapidly during the first few iterations and remains at a very low level for the rest of the simulation time (CCs and pro-tumoral cells have been mostly eliminated). This causes the total Hamiltonian of the TME (shown in green) to be aligned with the anti-tumoral Hamiltonian. In other words, a TME-Hamiltonian exhibiting a negative sign during the simulation shows that the system is in the Elimination phase, while the magnitude of the Hamiltonian indicates the energy of the agents’ interactions.

In the Equilibrium phase (Fig. 8), relative strengths of 0.3 versus 0.1 produce initially about 80 anti-tumoral cells and 60 pro-tumoral cells. In this case, a moderate increase in the population of pro-tumoral cells can be observed during the beginning of the simulations. As time passes this growth is counteracted by the anti-tumoral agents and the population of pro-tumoral cells remains quite stable for the remaining time of the simulations. Interestingly, as the number of anti-tumoral cells decreases substantially towards the end of the simulations, the number of pro-tumoral cells begins to show a moderate increase and the dispersion of the data grows larger. We believe that this behavior corresponds well with what is expected during the Equilibrium phase: in general and on average, the CCs are kept in check, but after balance has been maintained for a while the confrontation must be resolved either by the elimination of the CCs, or by their proliferation (leading to the Escape phase). Since the model contains stochastic elements, the resolution sometimes favors one side and sometimes favors the opposite, which is shown by the dispersion of the data. Values around zero of the corresponding TME-Hamiltonian reflect the balance between the opposing strengths from the beginning to the end of the simulations, even when the energy of the anti-tumoral and the pro-tumoral agents is not zero at the starting point.

In the Escape phase (Fig. 9), strength levels of 0.4 versus 0.6 generate initial populations of approximately 80 anti-tumoral cells and 100 pro-tumoral cells, respectively. In this case the pro-tumoral cells rapidly increase their numbers and surpass those of the anti-tumoral cells, although only up to a certain point (reaching about 110 cells on average). Then the population of pro-tumoral cells declines at varying rates up to the middle point of the simulations and grows again at a steady pace for the remaining time. Quite interestingly, the dispersion of the pro-tumoral data is significantly larger through the Escape phase than through the other two phases of CI. This is not caused by a larger spread in the sampling of the initial parameters (since the spread hyperparameter is the same as for the other phases), but is a result of the stochastic interaction in the presence of a larger population of pro-tumoral cells. The data of the TME-Hamiltonian also exhibits a larger dispersion, but a clear indication that the pro-tumor agents dominate the interaction in the TME is the fact that the Hamiltonian always shows a positive sign. In other words, a TME-Hamiltonian with a positive sign throughout the simulation indicates that the system is in the Escape phase. Also, in this phase the TME-Hamiltonian is always aligned with the pro-tumoral Hamiltonian and the magnitude of the energy is similar to that in the Elimination phase, but of opposite sign.

The results in Figs. 10, 11 and 12 show that the characteristic system behavior discussed above and corresponding to each of the CI phases can be obtained for different sets of hyperparameters: in these experiments the results were obtained using mean values for the anti-tumoral and the pro-tumoral agents with a constant difference between them (according to the equation reported in the caption of each of the referred figures). In other words, the hyperparameters’ values for which our model generates a certain phase of CI are neither unique nor arbitrary, but correspond to relationships between the values chosen for the anti-tumoral and the pro-tumoral agents that can be logically inferred from the desired system behavior. These results also demonstrate that the model is robust (not over-sensitive) to the hyperparameters, since small changes in these do not lead to qualitatively different responses, but to proportional changes within a qualitatively similar behavior.

The probability densities in Fig. 13 provide a visualization of the energy distribution for the latter half of the simulations. Rather than showing the system’s dynamics through time, these plots summarize the energy in the system as it reaches a stable state. By including five different hyperparameter settings (corresponding to those used to produce the plots in Figs. 10, 11 and 12), each of them simulated 50 times, we obtain enough data to generate conclusions with confidence. All the densities show that the highest frequency of occurrence is located around zero, except for the pro-tumoral energy (and consequently the TME energy) in the Escape phase. These results demonstrate how the system’s behavior is clearly captured by the energy densities as quantitative information that can be subjected to statistical analysis. For instance, examining the difference between the mean and the median values (Table 5) we can determine whether the corresponding distribution is left-skewed or right-skewed. A pro-tumoral median value of zero and a left-skewed TME distribution describe a system in the Elimination phase, while a TME median value of zero with a right-skewed distribution characterizes a system in the Equilibrium phase. A bimodal right-skewed TME distribution with a large median, signals that the system is in the Escape phase.

The densities of the TME in Fig. 13 are presented in a combined form in Fig. 14. This can be interpreted as snapshots of a system’s energy as it goes through the three phases of CI. In Elimination, the median pro-tumoral energy is zero. A system in Equilibrium has a median TME energy of zero but with a right-skewed distribution that indicates that CCs are kept in check but they are not being completely eliminated. A system in the Escape phase exhibits a lot of energy, both anti-tumoral and pro-tumoral; this can be related to the presence of more CCs compared to the other CI phases and more interaction between CCs and IS elements. The much wider spread of the distribution shows that systems in the Escape phase behave much more heterogeneously than in the other phases of CI. The bimodality of the distribution in the Escape phase indicates the presence of anti-tumoral and pro-tumoral elements in the system, while the higher frequency and larger values of the positive mode identify the prevalence of the CCs in the system.

Fig. 14
figure 14

Energy probability density estimates in different phase of CI


By conceiving the cellular-level interaction between cancer and the immune response as a confrontation between two adversaries (each of them with multiple agents) whose actions are governed by probabilistic rules, we developed a computational model of the biological system in the TME and framed this model under the more general Ising model from statistical mechanics. This allowed us to analyze the system’s status through time by means of its Hamiltonian. Through simulation of our proposed model in a multi-agent environment, we were able to verify its operation and generate the different phases involved in Cancer Immunoediting. These results show the ability of this model to capture complex stochastic behavior in a manner that is useful to understand the system and, eventually, control its outcome. Significantly scaling up the populations of cancer cells and IS elements is a future challenge that implies very time-consuming simulations, but we foresee that the results and conclusions would remain qualitatively similar to those reported herein.

The proposed quantitative analysis method integrates variables and parameters that model some of the biochemical elements in the TME; the system’s Hamiltonian captures these in an equation that expresses how the interaction among such simple cells emerges as the Cancer-IS complex behavior. The biological system in the TME is well described by our computational model, which can be tuned to general (Elimination, Equilibrium or Escape) or specific conditions of a unique organism, given the appropriate data. The cellular interaction in CI comprises stochastic processes that are accounted for in our current model by means of flexible probabilistic rules and driven in this study by Gompertzian growth. Nevertheless, the methodological proposal is not restricted to these choices; thus one can incorporate the formalisms that best fit the laboratory data available. It is worth mentioning that to date we have observed a scarcity of real data that could be employed in the validation of CI models; this unavailability of data is attributed to the complexity of the CI process and the numerous associated obstacles for its observation in the laboratory.

The method presented herein opens alternatives for exploring different aspects of the Cancer Immunoediting process; for instance, the variables and parameters corresponding to the biochemical elements in anti-cancer therapies may be directly integrated through the Ising-model Hamiltonian formalism, such that complementary extensions of the proposed model may capture the system’s dynamics that quantify the effect of diverse (chemo-, radio- or immuno-) therapies. In future work we will pursue these ideas by combining our findings with important contributions in related studies, particularly those related to DE-based and hybrid modeling of therapeutic elements. Our ultimate goal is the development of a comprehensive energy-based model of cancer immunology.

Availability of data and materials

The datasets generated and/or analyzed during the current study are not publicly available due to limitations imposed by our funding agency, but are available from the corresponding author on reasonable request.


  1. In computer science, an energy-based model associates a scalar probability (called energy) to each configuration of observed and latent variables; the name is imported from statistical physics to machine learning, without the literal meaning.

  2. This expression is found in Additional file 1: S1 Text of the cited work.



Agent-based model


Adaptive immune system


Cancer cells


Cancer immunoediting


Innate immune system


Immune system






Natural Killers


Ordinary differential equations


Partial differential equations


Stochastic differential equations


Tumor micro-environment


  1. Dunn GP, Bruce AT, Ikeda H, Old LJ, Schreiber RD. Cancer immunoediting: from immunosurveillance to tumor escape. Nat Immunol. 2002;3(11):991–8.

    Article  CAS  PubMed  Google Scholar 

  2. Dunn GP, Koebel CM, Schreiber RD. Interferons, immunity and cancer immunoediting. Nat Rev Immunol. 2006;6(11):836–48.

    Article  CAS  PubMed  Google Scholar 

  3. Narendra BL, Reddy KE, Shantikumar S, Ramakrishna S. Immune system: a double-edged sword in cancer. Inflamm Res. 2013;62(9):823–34.

    Article  CAS  Google Scholar 

  4. Mittal D, Gubin MM, Schreiber RD, Smyth MJ. New insights into cancer immunoediting and its three component phases-elimination, equilibrium and escape. Curr Opin Immunol. 2014;27:16–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Hanahan D, Weinberg R. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74.

    Article  CAS  PubMed  Google Scholar 

  6. Norton K-A, Gong C, Jamalian S, Popel A. Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes. 2019;7(1):37.

    Article  CAS  PubMed  Google Scholar 

  7. Mahlbacher GE, Reihmer KC, Frieboes HB. Mathematical modeling of tumor-immune cell interactions. J Theor Biol. 2019;469:47–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Suhail Y, Cain MP, Vanaja K, Kurywchak PA, Levchenko A, Kalluri R, et al. Systems biology of cancer metastasis. Cell Syst. 2019;9(2):109–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Esquivel-Velázquez M, Ostoa-Saloma P, Palacios-Arreola MI, Nava-Castro KE, Castro JI, Morales-Montor J. The role of cytokines in breast cancer development and progression. J Interferon Cytokine Res. 2015;35(1):1–16.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Rhodes A, Hillen T. A mathematical model for the immune-mediated theory of metastasis. J Theor Biol. 2019;482:109999.

    Article  CAS  PubMed  Google Scholar 

  11. Barriga V, Kuol N, Nurgali K, Apostolopoulos V. The complex interaction between the tumor micro-environment and immune checkpoints in breast cancer. Cancers. 2019;11(8):1205.

    Article  CAS  PubMed Central  Google Scholar 

  12. Lim HX, Kim TS, Poh CL. Understanding the differentiation, expansion, recruitment and suppressive activities of myeloid-derived suppressor cells in cancers. Int J Mol Sci. 2020;21(10):3599.

    Article  CAS  PubMed Central  Google Scholar 

  13. Kumar V, Patel S, Tcyganov E, Gabrilovich DI. The nature of myeloid-derived suppressor cells in the tumor microenvironment. Trends Immunol. 2016;37(3):208–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Szeto GL, Finley SD. Integrative approaches to cancer immunotherapy. Trends Cancer. 2019;5(7):400–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Yang L, Pang Y, Moses HL. TGF-β and immune cells: an important regulatory axis in the tumor microenvironment and progression. Trends Immunol. 2010;31(6):220–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Dahmani A, Delisle J-S. TGF-\(\beta\) in T cell biology: implications for cancer immunotherapy. Cancers. 2018;10(6):194.

    Article  PubMed Central  CAS  Google Scholar 

  17. Peinado H, Zhang H, Matei IR, Costa-Silva B, Hoshino A, Rodrigues G, Psaila B, Kaplan RN, Bromberg JF, Kang Y, Bissell MJ, Cox TR, Giaccia AJ, Erler JT, Hiratsuka S, Ghajar CM, Lyden D. Pre-metastatic niches: organ-specific homes for metastases. Nat Rev Cancer. 2017;17(5):302–17.

    Article  CAS  PubMed  Google Scholar 

  18. Pulido C, Vendrell I, Ferreira AR, Casimiro S, Mansinho A, Alho I, Costa L. Bone metastasis risk factors in breast cancer. Ecancermedicalscience. 2017;11:715.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Gabrilovich DI, Nagaraj S. Myeloid-derived suppressor cells as regulators of the immune system. Nat Rev Immunol. 2009;9(3):162–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. McCoy BM, Wu TT. The two-dimensional ising model. Cambridge, MA: Harvard University Press; 2013.

    Google Scholar 

  21. Barradas-Bautista D, Alvarado-Mentado M, Agostino M, Cocho G. Cancer growth and metastasis as a metaphor of Go gaming: an ising model approach. PLoS ONE. 2018;13(5):0195654.

    Article  CAS  Google Scholar 

  22. Rojas-Dominguez A, Barradas-Baustista D, Alvarado M. Modeling the game of go by ising Hamiltonian, deep belief networks and common fate graphs. IEEE Access. 2019;7:120117–27.

    Article  Google Scholar 

  23. Lander A. Pattern, growth, and control. Cell. 2011;144(6):955–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Werner HM, Mills GB, Ram PT. Cancer systems biology: a peek into the future of patient care? Nat Rev Clin Oncol. 2014;11(3):167.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Pourhasanzade F, Sabzpoushan SH, Alizadeh AM, Esmati E. An agent-based model of avascular tumor growth: immune response tendency to prevent cancer development. Simulation. 2017;93(8):641–57.

    Article  Google Scholar 

  26. Jinnah A, Zacks B, Gwam C, Kerr B. Emerging and established models of bone metastasis. Cancers. 2018;10(6):176.

    Article  PubMed Central  CAS  Google Scholar 

  27. Denton AE, Roberts EW, Fearon DT. Stromal cells in the tumor microenvironment. In: Owens BMJ, Lakins MA, editors. Stromal immunology, vol. 1060. Cham: Springer; 2018. p. 99–114.

    Chapter  Google Scholar 

  28. Shiozawa Y, Eber MR, Berry JE, Taichman RS. Bone marrow as a metastatic niche for disseminated tumor cells from solid tumors. BoneKEy Rep. 2015;4:689.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Méndez-García LA, Nava-Castro KE, Ochoa-Mercado TL, Palacios-Arreola MI, Ruiz-Manzano RA, Segovia-Mendoza M, Solleiro-Villavicencio H, Càzarez-Martínez C, Morales-Montor J. Breast cancer metastasis: are cytokines important players during its development and progression? J Interferon Cytokine Res. 2019;39(1):39–55.

    Article  PubMed  CAS  Google Scholar 

  30. Altrock PM, Liu LL, Michor F. The mathematics of cancer: integrating quantitative models. Nat Rev Cancer. 2015;15(12):730–45.

    Article  CAS  PubMed  Google Scholar 

  31. Gong C, Milberg O, Wang B, Vicini P, Narwal R, Roskos L, Popel AS. A computational multiscale agent-based model for simulating spatio-temporal tumour immune response to PD1 and PDL1 inhibition. J R Soc Interface. 2017;14(134):20170320.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Newton PK, Mason J, Bethel K, Bazhenova L, Nieva J, Norton L, Kuhn P. Spreaders and sponges define metastasis in lung cancer: a Markov chain Monte Carlo mathematical model. Can Res. 2013;73(9):2760–9.

    Article  CAS  Google Scholar 

  33. Cohen E, Gao H, Anfossi S, Mego M, Reddy N, Debeb B, et al. Inflammation mediated metastasis: immune induced epithelial-to-mesenchymal transition in inflammatory breast cancer cells. PLoS ONE. 2015;10:e0132710.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Figueredo GP, Aickelin U. Defining a simulation strategy for cancer immunocompetence. In: International conference on artificial immune systems. Springer; 2010. p. 4–17.

  35. Novozhilov AS, Berezovskaya FS, Koonin EV, Karev GP. Mathematical modeling of tumor therapy with oncolytic viruses: regimes with complete tumor elimination within the framework of deterministic models. Biol Direct. 2006;1(1):1–18.

    Article  CAS  Google Scholar 

  36. Chiacchio F, Pennisi M, Russo G, Motta S, Pappalardo F. Agent-based modeling of the immune system: NetLogo, a promising framework. Biomed Res Int. 2014;2014:1–6.

    Article  Google Scholar 

  37. Wilson S, Levy D. A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy. Bull Math Biol. 2012;74:1485–500.

    Article  CAS  PubMed  Google Scholar 

  38. Shahriyari L. A new hypothesis: some metastases are the result of inflammatory processes by adapted cells, especially adapted immune cells at sites of inflammation. F1000Research. 2016;5:175.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Blomberg OS, Spagnuolo L, de Visser KE. Immune regulation of metastasis: mechanistic insights and therapeutic opportunities. Dis Models Mech. 2018;11(10):036236.

    Article  CAS  Google Scholar 

  40. Newton PK, Mason J, Bethel K, Bazhenova LA, Nieva J, Kuhn P. A stochastic Markov chain model to describe lung cancer growth and metastasis. PLoS ONE. 2012;7(4):34637.

    Article  CAS  Google Scholar 

  41. Newton PK, Mason J, Hurt B, Bethel K, Bazhenova L, Nieva J, Kuhn P. Entropy, complexity and Markov diagrams for random walk cancer models. Sci Rep. 2015;4(1):7558.

    Article  CAS  Google Scholar 

  42. Newton PK, Mason J, Venkatappa N, Jochelson MS, Hurt B, Nieva J, Comen E, Norton L, Kuhn P. Spatiotemporal progression of metastatic breast cancer: a Markov chain model highlighting the role of early metastatic sites. NPJ Breast Cancer. 2015;1(1):15018.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Wells DK, Chuang Y, Knapp LM, Brockmann D, Kath WL, Leonard JN. Spatial and functional heterogeneities shape collective behavior of tumor-immune networks. PLoS Comput Biol. 2015;11(4):1004181.

    Article  CAS  Google Scholar 

  44. Ku-Carrillo RA, Delgadillo SE, Chen-Charpentier B. A mathematical model for the effect of obesity on cancer growth and on the immune system response. Appl Math Model. 2016;40(7–8):4908–20.

    Article  Google Scholar 

  45. Núñez-López M, Hernández-López E, Delgado J. Stochastic simulation on a minimal model of cancer immunoediting theory. Int J Bifurc Chaos. 2021;31(06):2150088.

    Article  Google Scholar 

  46. Figueredo GP, Siebers P-O, Aickelin U. Investigating mathematical models of immuno-interactions with early-stage cancer under an agent-based modelling perspective. BMC Bioinform. 2013;14(6):1–20.

    Google Scholar 

  47. Figueredo GP, Siebers P-O, Owen MR, Reps J, Aickelin U. Comparing stochastic differential equations and agent-based modelling and simulation for early-stage cancer. PLoS ONE. 2014;9(4):95150.

    Article  CAS  Google Scholar 

  48. Sajid N, Convertino L, Friston K. Cancer niches and their kikuchi free energy. Entropy. 2021;23(5):609.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Torquato S. Toward an ising model of cancer and beyond. Phys Biol. 2011;8(1):015017.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Wilkie KP, Hahnfeldt P. Modeling the dichotomy of the immune response to cancer: cytotoxic effects and tumor-promoting inflammation. Bull Math Biol. 2017;79(6):1426–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Cleveland C, Liao D, Austin R. Physics of cancer propagation: a game theory perspective. AIP Adv. 2012;2(1):011202.

    Article  PubMed Central  Google Scholar 

  52. Mahlbacher G, Curtis LT, Lowengrub J, Frieboes HB. Mathematical modeling of tumor-associated macrophage interactions with the cancer microenvironment. J Immunother Cancer. 2018;6(1):10.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Lewkiewicz S, Chuang Y-L, Chou T. A mathematical model of the effects of aging on naive T cell populations and diversity. Bull Math Biol. 2019;81(7):2783–817.

    Article  PubMed  Google Scholar 

  54. Weber M, Buceta J. The cellular ising model: a framework for phase transitions in multicellular environments. J R Soc Interface. 2016;13(119):20151092.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Wilensky U. NetLogo. Evanston, IL: Center for Connected Learning and Computer-Based Modeling, Northwestern University; 1999. Accessed 15 Nov 2020.

  56. Peter DH. Kernel estimation of a distribution function. Commun Stat-Theory Methods. 1985;14(3):605–20.

    Article  Google Scholar 

Download references


The authors acknowledge Mariana Arroyo Duarte, M.D., pathologist with the National Mexican Institute of Social Security (IMSS), for her valuable advice regarding medical aspects discussed in this paper. Daniela Flores Silva was involved in the development and implementation of an earlier version of the model presented herein. Sergio Alcalá helped in the production of simulation plots for Additional file 1. The authors acknowledge their contribution.


This work was supported by the National Council of Science and Technology of Mexico (CONACYT) through grants: A1-S-20037 (M. Alvarado, lead researcher) and CÁTEDRAS-2598 (A. Rojas-Domínguez).

Author information

Authors and Affiliations



ARD participated in the design of the study and algorithms, and in the implementation of the model, carried out computational experiments, analyzed experimental results, and wrote the paper with contributions from the other authors. RAD participated in the formal characterization of the interactions in cancer immunoediting and implemented initial versions of the model. FRV implemented the model and carried out computational experiments. MA conceived of the study and the essential elements of the model, participated in the design of the algorithms, analyzed the results, participated all along in the preparation of the manuscript, and procured the publication fees. All authors read and approved the final manuscript.

Authors’ Information

Alfonso Rojas-Domínguez holds M.Sc. and Ph.D. degrees from the University of Liverpool, UK. Since 2014, he has been a Research Fellow with the National Council of Science and Technology of Mexico and is a member of the National Researchers System in that country. His main interests are: machine learning, deep learning, computational intelligence and computational optimization, with biomedical applications.

Renato Arroyo-Duarte received an M. Sc. in physics from the University of Guadalajara and is currently a Ph.D. candidate at the same institution. His doctoral work is focused on the implementation of statistical physics methodologies in biological systems. He has also worked on information theory with applications in the study of genomic sequences.

Fernando Rincón-Vieyra is a Bachelor of Engineering in Computer Systems from IPN in Mexico; his scientific interest are systems medicine, machine learning and data science for healthcare. He is currently pursuing posgraduate studies.

Matías Alvarado holds a Ph.D. degree in mathematics from the Technical University of Catalonia and is a leader researcher in the Department of Computer Science, CINVESTAV-IPN. His ongoing research focuses on the formal modeling and algorithmic simulation of cancer metastasis and the immune system response. Previously he has studied selection strategies for sport team games using Nash equilibrium, vehicle’s speed adaptation through machine learning, knowledge management, and epistemic logic.

Corresponding author

Correspondence to Matías Alvarado-Mentado.

Ethics declarations

Ethics approval and consent to participate

Not Applicable.

Consent for publication

Not Applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Alfonso Rojas-Domínguez: CONACYT Research Fellow

Supplementary Information

Additional file 1.

Sensitivity analysis and analysis of model elements.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rojas-Domínguez, A., Arroyo-Duarte, R., Rincón-Vieyra, F. et al. Modeling cancer immunoediting in tumor microenvironment with system characterization through the ising-model Hamiltonian. BMC Bioinformatics 23, 200 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: