Agent-based modeling for the tumor microenvironment (TME)

1.   Introduction

Cancer is a disease that involves uncontrolled abnormal (tumor) cell growth. Due to this uncontrolled growth, the tumor cell density increases, which leads to a cell mass known as a tumor. The disease is usually incurable if the tumor cells spread to the surrounding tissues (invasion) or other organs (metastasis) (see Figure 1).

The tumor microenvironment (TME) is the environment that surrounds a tumor and consists of different types of cells such as cancer associated fibroblasts (CAF), cancer associated macrophages (CAM), T-cells, B-cells, pericytes, tumor cells, and the extracellular matrix (ECM) ^[1,2] (see Figure 2). Various components of the TME promote tumor progression, invasion, and metastasis ^{[3,4,5,6,7,8]}. Understanding how tumor cells adapt to their environment is one of the biggest challenges in developing successful cancer treatments. This research study seeks to understand the tumor cell behavior in the TME by investigating cell-cell and cell-ECM interactions through an agent-based model.

The physical barrier for all cells in the TME is the ECM. It contains angiogenic factors and chemokines, which give each tissue its tensile, comprehensive strength and elasticity ^[9]. Different proteins in the ECM make a surrounding membrane and an interstitial matrix ^{[10,11,12,13]}. Cell-ECM adhesion is most important for cell movement through the TME. ECM stiffness interrupts drug delivery and promotes metastasis ^[10,14]. When a tumor cell attaches to an ECM component, the tumor cells secrete collagenase proteases called matrix metalloproteinases (MMPs), which break down the ECM component ^{[15,16,17,18]}. The ECM stores pro-angiogenic growth factors that support the formation of new blood vessels close to the tumor cells ^[19,20].

Restricted nutrient or oxygen supply in the TME promotes resource competition between cells. Tumor cells follow different strategies such as scavenging ECM proteins via integrins, receptor-mediated albumin uptake, catabolism, macropinocytosis consumption of multiple components of the TME, and engulfment and degradation of entire live cells via entosis to increase their nutrient supply ^[21,22]. The metabolic reprogramming mechanism called the Warburg effect increases the energy of tumor cells more than healthy cells ^{[23,24,25,26,27]}.

Cell competition can occur through the competition for nutrients or growth factors and the competition for space through mechanisms such as cell fitness sensing ^{[28,29,30,31,32]}. The fitness level of a cell determines the proliferation and apoptosis of a cell. Super competition is the process of eliminating healthy wild-type cells with tumor cells. Myelocytomatosis (MYC) protein and Yes-associated proteins (YAP) are two proteins that turn cells into so-called super competitors ^[33]. In MYC-mediated cell competition, the losers with a low level of MYC are killed by the winners with a high level of MYC ^[34]. Tumor cells with metabolic winner phenotypes compete with immune cells and escape from the TME ^[35]. Hence, cell competition promotes tumor progression by killing less fit cells in the TME.

An epithelial to mesenchymal transition (EMT) occurs due to the disruption of cell-cell adhesion and cellular polarity, remodeling of the cytoskeleton, and changes in the cell-matrix adhesion ^[3]. In the deleterious adaptation of cells, cell-cell adhesion can influence how cells either escape or remain restricted to certain domains ^{[36,37,38,39]}.

Clinical experiments are time-consuming, expensive, and have issues capturing the system dynamics. However, some advantages of mathematical modeling are the ability to identify critical factors and the most sensitive parameters of a system, the ability to modify the model according to the objectives of the study, and the ability to predict the behavior of a system. Additionally, mathematical models can capture the dynamics of systems in time and space, thus enabling us to make sound analyses and predictions on which new treatment strategies can be developed.

Different mathematical models have been developed in the literature to analyze tumor growth, cancer metabolism, vasculature, tumor immunity, tumor microenvironment and heterogeneity, invasion, metastasis, and cancer treatments ^{[40,41,42,43]}. Basically, three modeling techniques have been used to understand cancer dynamics: discrete (agent-based), continuum, and hybrid. Agent-based models track and update individual cells according to a set of biological rules ^[44,45]. Through an agent-based model, one can link tumor heterogeneity to a cellular response. Furthermore, we can study how molecular signaling, cellular metabolism, mutation-induced phenotypic changes, angiogenesis, metastasis, immune response, and cancer treatments affect cancer progression ^[46,47]. Continuum models consider the tumor tissue as a continuous medium, and differential equations are used for modeling ^{[48,49,50,51,52,53]}. Diffusion models for tumors have been developed under different mechanisms such as negative feedback ^[54], discontinuous switches ^[55,56], non-uniformly distributed mitotic inhibitors ^[57,58], and through the effects of a necrotic core ^[59]. Compared to discrete models, it is easier to estimate parameters for continuum models based on laboratory experiments. On the other hand, some assumptions used in continuum modeling can be too simplistic. A hybrid model combines the benefits of discrete and continuous modeling and techniques and describes chemical reactions and tissue landscapes in a single model ^[60,61]. Hybrid or continuum-discrete models can bridge the gap between the cellular scale and the tumor scale. In the hybrid approach, tumor tissue is modeled, and contains discrete and continuum elements. Oxygen, nutrient, drug, growth factors, and certain tissue features are described as continuum fields, and cells are described as discrete elements.

Numerous agent-based models have been developed for cancer modeling, thereby incorporating a wide range of features such as proliferation, cytotoxic behaviors, mechanical pressure, cell death, cell movement, treatments, mutations, angiogenic factors, immune system response, and energy availability ^{[18,46,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86]}. For instance, Lettort et al. introduced an open-source package called PhysiBoss ^[87,88], which integrated intracellular signaling via Boolean modeling (MaBoss) ^[89] with multicellular behavior through agent-based modelling (PhysiCell) ^[90]. This approach is particularly advantageous to study the responses of heterogeneous cell populations to treatments, mutation effects, and different modes of invasion or isomorphic morphogenesis events ^[91]. However, this model does not account for ECM interactions within the tumor microenvironment (TME).

Previously developed agent-based models often only consider a few features of the TME. In order to understand the role of the TME in a greater depth, we focus on the eight hallmarks in our model suggested by Hanahan and colleagues ^[92,93,94]. The key hallmarks are given as follows:

● Sustaining proliferative signaling: Tumor cells can continuously signal themselves to grow, thus bypassing the normal regulatory mechanisms.

● Evading growth suppressors: Tumor cells can evade the growth-suppressing signals that would normally inhibit uncontrolled cell division.

● Enabling replicative immortality: Tumor cells maintain their telomeres, thus allowing them to divide indefinitely, unlike healthy cells that have a limited number of divisions.

● Resisting cell death: Tumor cells develop mechanisms to avoid programmed cell death, enabling their survival and proliferation.

● Deregulating cellular energy: Tumor cells reprogram their energy production to fuel rapid growth and adapt to the tumor microenvironment.

● Avoiding immune destruction: Tumor cells develop mechanisms to evade immune detection and destruction.

● Inducing angiogenesis: Tumors can stimulate the growth of new blood vessels to supply nutrients and oxygen, which are critical for their continued growth.

● Invasion and metastasis: Tumor cells can invade surrounding tissues and spread to distant organs, thus forming secondary tumors.

Figure 3 shows features of the model related to each cancer hallmark. In this study, we develop an agent-based model to study interactions between the tumor cells and the TME according to a set of biological rules. We emphasise that our model does not account for the formation of blood vessels. We acknowledge that the existence of blood vessels can play an important role in the tumor spread to the boundary of an organ and that the location of the tumor in the body is important; however, since we do not account for the blood vessel location, this important issue is beyond the scope of our model. Our aim is to keep the model simple but insightful. Furthermore, we observe the time that the first tumor cell exits the TME under various environmental conditions. This first cell exit time measures how long a tumor takes to spread to the surrounding tissues and other organs, and this information can help in the development of successful treatment strategies, as cancers are life-threatening when they invade or metastasize.

In Section 2, we describe the agent-based model that we developed to study the tumor cell behavior at the cellular level. Simulation results of this agent-based model are described in Section 3. Section 4 provides a discussion, and conclusions are provided in Section 5.

2.   Agent-based model description

This model represents a two-dimensional tissue section that consists of healthy cells, tumor cells, and ECM proteins. The domain is a $100\times 100$ grid with two regions (inner and outer regions). The outer region of the domain represents the invasive margin of the tissue. We assume that the domain is initially a healthy tissue (see Figure 4).

This tissue section develops into a tumor due to a series of mutations of a healthy cell. We assume that healthy cells are only located in the inner region and their movements are negligible as adhesion molecules keep the healthy cells in their intended locations ^[37,95]. We randomly locate ECM proteins in the domain and assume that they cannot move. Furthermore, tumor cells cannot occupy a grid position of an ECM protein. The boundary of the domain has no gaps, is said to be fully fenced at time zero, and represents the surrounding membrane of the tissue. The maximum population in the environment is 10,000, as we have a total of $10,000$ grid positions in the domain.

The input parameters and the calculated output variables for the agent-based model are listed in Table 1.

This agent-based model observes the effect of different factors of the TME on tumor progression, plasticity, and invasion. The output variables $N(t)$ and $A(t)$ (healthy and tumor cell densities, respectively, at time $t$) measure tumor progression under different circumstances. We consider two output variables, $\bar{S}(t)$ and $\bar{J}(t)$, to represent the effect of the tumor microenvironmental factors on the tumor cell plasticity. To measure the effect of different factors of the TME for tumor invasion, we use the two output variables, $E(t)$ and $f(t)$ - see Table 1.

Figure 5 represents the flow diagram of the model and we describe the set of rules associated with the agent-based model. We select a tumor cell in the domain and check the number of divisions of the cell. When a tumor cell has completed the maximum number of cell divisions of a healthy cell ($K$), it has a chance of dying. If a cell dies, then the grid position that the cell occupied becomes empty. If the cell's age is greater than the division age, then that cell has a chance to proliferate. When a cell proliferates, the daughter cell is placed in the first empty grid position in a clockwise direction starting from the North (using a Moore neighborhood). If all neighboring grid positions are occupied, then the mother cell competes with the weakest cell in the neighborhood determined by the competition rate, and that position is used by the daughter cell (see Figure 6). At any age, a cell can move to a new position depending on the number of cells and/or obstacles in the neighborhood and the cell stickiness. We assume that a tumor cell can move North, East, South or West (movement is restrained to a Von Neumann neighborhood). A detailed description of the model parameters, output variables, and rules is given in the Appendix.

The results of this agent-based model provide a good understanding of the individual tumor cell behavior under different conditions of the TME. Through this study, we can attempt to identify the most influential tumor microenvironmental factors that affect the tumor progression, cell plasticity, and invasion.

3.   Simulation results of the agent-based model

To study tumor cell behavior in the TME, we perform numerical simulations of the agent-based model using the MATLAB programming language. The relevant program codes can be accessed through https://github.com/Hasi1988/abm_for_TME. We observe the changes of the cell densities, average cell stickiness, average jump radii and fence densities with time. Furthermore, we observe the total number of tumor cells that exit from the domain. In this study, we also observe the effect of three factors on the TME: the initial healthy cell densities ($N(0)$), ECM protein densities ($R$) and ECM breakdown probabilities ($P_b$). We explore tumor progression, tumor plasticity, and tumor invasion under changes in each of these three factors. For each scenario, we perform $500$ simulations. To investigate the effect of a factor, we change just that factor and keep the other parameter values, as in Table 2. For example, to investigate the effect of $N(0)$, we only consider $N(0) = 0.243, 0.324,$ and $0.405,$ while keeping $R = 0.1$ and $P_b = 0.5$. Other scenarios consider values of $0.25$ and $0.75$ as the ECM breakdown probabilities, and values of $0.05$ and $0.2$ as the obstacle densities.

We assume that the healthy cells are only located in the inner region that consists of $81\times81$ grid positions. We consider that the degradation of the ECM by a tumor cell can be equally likely to either occur or not to occur by assuming the ECM breakdown probability as $0.5$. Since the protein density differs from tissue to tissue ^[96], our model assumes that the ECM protein density is $0.1$. Due to adhesion molecules, the healthy cells have a high stickiness and a low motility ^[37]. Hence, we assume that the cell stickiness and jump radius of healthy cells are $0.9$ and $1$, respectively. In 1965, Hayflick suggested that healthy cells can only divide a limited number of times by observing the division times of human diploid fibroblast cells ^[97]. However, damage signals such as DNA-damaging agents and oxidative stress can affect the early cessation of the division of cells ^[98]. Since the TME is a hypoxic environment, we assume that a healthy cell can divide a maximum of $10$ times. Our model does not account for the formation of blood vessels, so whether angiogenesis has occurred is not relevant for our model. We wish to keep our model as simple as possible. We consider a dimensionless time from $0$ to $24$ by dividing this into $1800$ sub-intervals. Hence, the time step length is $0.0133 \:\:(24/1800)$. We assume that the age of each cell increases by $1$ at each time step and the division age of a cell is $20$. Since tumor cells can proliferate uncontrollably, we assume a high probability rate ($0.6$) for cell division. Furthermore, we assume a significantly small probability of death ($0.001$) due to the resistance to death of tumor cells ^[93]. Since the complexity of the TME and the condition of the environment can differ from patient to patient, the parameter values can vary according to the condition of the TME.

Figure 7 represents the progression of a tumor in healthy tissue. With uncontrolled growth, the tumor cell density increases by creating a crowded environment. Due to the lack of nutrients and space, tumor cells compete with the surrounding cells. Since the energy of healthy cells decreases with the number of divisions ^[99,100], most healthy cells are not strong enough to compete with super-competitive tumor cells. Hence, many healthy cells in the environment die and provide space for tumor cells. Moreover, tumor cells regulate the escape mechanisms to move from this crowded, stressful environment. By secreting MMPs, aggressive tumor cells degrade the ECM components in the boundary. At the next step, malignant tumor cells escape through the boundary and invade the surrounding tissues.

Subsections 2.1–2.3 describe the scenario outcomes when varying the initial healthy cell density, the ECM protein density, and the ECM breakdown probability for the tumor progression, plasticity, and invasion, respectively.

3.1.   Tumor progression under the effect of $N(0)$, $R$ and $P_b$

When a tumor develops, the healthy cell density of the environment decreases and the tumor cell density increases. We observe the healthy and tumor cell densities for each scenario to investigate the tumor progression.

We can observe from the scenarios exploring the tumor progression (see Figures 8, 9 and Table 3) through the following aspects:

● A larger number of initial healthy cells slightly slows the increase in tumor cell density and results in a slightly lower tumor cell density at equilibrium.

● The density of the healthy cells at equilibrium is reduced (due to tumor growth), but is closely related to the initial number of healthy cells.

● The density of obstacles (or proteins) in the ECM impacts the growth of the tumor. At equilibrium, the tumor cell density is significantly less for an obstacle density of 0.2 compared to a density of 0.05.

● Accordingly, more healthy cells survive at equilibrium if the ECM has a higher density of obstacles.

● It takes slightly longer for the tumor cell density to reach equilibrium for a larger initial healthy cell density or for a larger density of obstacles.

● Tumor cells compete for space and nutrients, and kill weak healthy cells. An environment with a low density of obstacles and a lower healthy cell density is advantageous for tumor cells. Conversely, the tumor progression is impacted in a crowded environment.

● Simulation results showed that the ECM breakdown probability does not significantly affect the tumor progression. The graphs were omitted as they largely overlapped, thus indicating minimal variation and a negligible impact on the tumor progression. This suggests that the ECM breakdown probability does not play a critical role in the progression of the tumor.

● Low variations at early times and after equilibrium show the homogeneous state of the environment. The high variation in the cell densities during the middle period suggests that the system is undergoing significant changes.

3.2.   Tumor cell plasticity under the effect of $N(0)$, $R$ and $P_b$

We investigate the ability of tumor cells to change their phenotypes according to their environment. In particular, we observe the change in the average cell stickiness and the average jump radius of the tumor cells. In these simulations, the initial stickiness value for a healthy cell is $0.9$. The stickiness of a tumor cell differs depending on the change in stickiness and the stickiness of the mother cell - see Table 1. We neglect the movements of the healthy cells, as they have a high stickiness and a low motility. The jump radius of a tumor cell changes according to the jump radius of the mother cell and the change in jump radius - again see Table 1. Figure 10 and Table 4 show the average stickiness of tumor cells. From Tables 3 and 4 and Figures 9 and 10, we can make the following comments:

● The average cell stickiness of the tumor cells reduces over time and is lower for increasing initial healthy cell densities. This stickiness value reduces to a lower value for higher values of obstacle densities.

● If the environment is crowded, then the tumor cells will lose their stickiness; hence, they will not be as likely to collectively move in a group, and this slows down the tumor progression.

● The average jump radius of the tumor cells increases more rapidly for lower initial healthy cell densities, although it approaches a similar equilibrium value regardless of $N(0)$.

● The average jump radius of the tumor cells increases more rapidly and approaches a higher equilibrium value when the obstacle density is lower.

● Tumor cells regulate the cell plasticity mechanisms according to the various conditions of the TME. In a crowded environment, the tumor cells can reduce their stickiness and increase their jump radius, which affects their ability to move. Hence, the tumor cells promote cell plasticity.

● The initial healthy cell density and the ECM protein density influence the tumor cell plasticity regulations. The results suggest that if there are many environmental barriers, then the tumor cells try to migrate as individual cells by reducing their stickiness at a high rate. Otherwise, they may escape the environment as a collection of cells. Hence, the invasive pattern of the cells may differ according to the tumor microenvironmental conditions.

● The tumor cells have high motility rates under low initial healthy cell densities and ECM protein densities. Furthermore, they may try to move as a collection of cells under those conditions. Hence, the high initial healthy cell densities and ECM protein densities in the TME slow down the motility rates of the tumor cells.

● The simulation results showed that the ECM breakdown probability did not significantly affect the tumor cell plasticity. The graphs were omitted as they are similar to results already presented. This indicates a minimal variation and a negligible impact on the tumor cell plasticity. This suggests that the ECM breakdown probability does not play a critical role in the tumor cell plasticity.

3.3.   Tumor cell invasion under the effect of $N(0)$, $R$ and $P_b$

We investigate the effects of the tumor microenvironmental factors that lead to an increased likelihood of tumor invasion. We measure this by determining the number of tumor cells that exit the tumor microenvironment, as well as the degradation of the fence density that provides greater chances of tumor invasion. See Figure 11 and average first cell exit times given in Table 5. We note the following:

● More tumor cells exit the TME when there is a reduced initial healthy cell density, as these conditions lead to faster tumor progression.

● More tumor cells exit the TME when the probability of a breakdown in the fence density of the ECM ($P_b$) increases.

● The density of the ECM proteins ($R$) impacts the tumor cell exit. More tumor cells exit the TME for an obstacle rate of $0.05$ compared to a rate of $0.2$.

● Accordingly, a lesser number of tumor cells exit for increased initial healthy cell densities, increased ECM proteins, and low ECM breakdown probabilities.

● A reduction of the fence density is slow for high initial healthy cell densities, low ECM breakdown probabilities, and high ECM protein densities.

● As the stickiness of a cell reduces to a lower value for the higher values of the initial healthy cell densities and obstacle densities (see Figure 10), many cells may individually exit the TME through existing gaps without degrading new fence components under these conditions.

● The first cell exit time is influenced by the density of the ECM obstacles/proteins; a lower density leads to an earlier first exit time.

● The breakdown probability does not have much influence on the first exit time, neither does the initial healthy cell density.

If a healthy tissue has many cells or ECM proteins, then the cell invasion speed is slow. Furthermore, the tumor cell invasion can slow if the ECM breakdown probability of the tumor cells is low. Moreover, the tumor cells start to exit the environment early for low initial healthy cell densities, high ECM breakdown probabilities, and low obstacle densities. We observe that the tumor progression is high under these tumor microenvironmental conditions. Hence, the tumor cells try to exit early, as the environment becomes quickly crowded under these conditions.

4.   Discussion

In this paper, we explored how various components of the tumor microenvironment promote tumor progression, invasion, and metastasis via a mathematical agent-based model. In this agent-based model, we simulated cell-cell and cell-ECM interactions according to a set of biological rules, underpinned by the eight hallmarks of cancer suggested by Hanahan and colleagues ^[92,93,94]. This allowed for a link between the tumor heterogeneity and the cellular response. Then, the output from the model suggested lines of enquiry into the system dynamics, thus augmenting clinical experiments that are often time-consuming and expensive. To get such insights into these dynamics (in time and space) is invaluable to predict possible new treatment strategies. We re-emphasize that we did not calibrating our model against experimental results. Rather, we adopted a standard approach of using detailed simulations to explore potential hypotheses.

Cancer is life-threatening if the tumor cells spread into the surrounding tissues and other organs. The results suggest that the ECM protein density, the initial healthy cell density, and the ECM degradation rate affect tumor cell invasion. Hence, by focusing on these three factors, we can explore the spread of tumor cell invasion. For example, the tumor cell invasion spread is slow if there are many ECM proteins and healthy cells in the domain. By using an agent-based model similar to the one described here, biologists, clinicians, and modelers can work together to test different hypotheses that can offer insights into the design of new therapies. In particular, this synergy between different groups allows the researchers to eliminate scenarios that are not appropriate in designing new cancer therapies.

Our agent-based model represents a two-dimensional tissue section that consists of healthy cells, tumor cells, and ECM proteins; the domain is a 100 $\times$ 100 grid that initially represents healthy tissue. The simulation started with one healthy cell mutating, which allowed a tumor to grow. By setting various parameters for the model (e.g., the density of the ECM proteins in the domain; the stickiness of the cells that affects the tumor cell movement; and the density of the fence around the boundary of the domain that affects the likelihood of a tumor cell exiting and spreading beyond the domain), the simulation tracked the density of healthy cells and tumor cells over time and thus provides a link between the parameters and tumor progression, invasion, and metastasis.

We observed that tumor progression was affected by the initial density of healthy cells in the domain and the density of the ECM proteins. We observed tumor cell plasticity in terms of a reduction in the tumor cell stickiness (that may be due to crowdedness of the environment) that influenced how the tumor cells escaped the domain and spread; additionally, this is associated with the probability of the ECM protein breakdown along the boundary of the domain.

Our model is a starting point to explore tumor growth, and a limitation is that it is only spatially two-dimensional and is based on a regular mesh. Consequently, it only handles a limited complexity, and there are still issues with how to associate the results from the simulations with clinical outcomes. Furthermore, there is clearly a need to carry out a sensitivity analysis; however, due to the number of parameters, this has not proved feasible in this present work. Additionally, it would be interesting to compare our lattice-based approach with off-lattice models, though this is beyond the scope of this present work.

5.   Conclusions

This paper explored the effects of the initial healthy cell density, the ECM protein density, and the ECM breakdown probability for tumor progression, cell plasticity, and invasion using an agent-based model. Although the model was based on a two-dimensional spatial domain, we are able to demonstrate links between the initial healthy cell densities, the ECM protein densities, and the cell stickiness with tumor progression, cell motility, and invasion. We can easily change the model parameters to gain insights into the tumor cell behavior. Additionally, it is possible to modify the model to study different treatment strategies, such as the effect of T-cells on cancer dynamics, which will be further investigated at a later date. Incorporating T-cell effects can contribute to additional knowledge and offer new insights into novel treatment protocols; therefore, the study in this paper can be considered as an initial step in this direction.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research work was supported by the Australian Research Council Future Fellowship FT180100698.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A description of input parameters and calculated output variables for the agent-based model

(1) Healthy Cell Density $(N(t))$ - We calculate $N(t)$ by dividing the current healthy cell population by total of $10,000$ grid positions.

(2) Tumor Cell Density $(A(t))$ - We calculate $A(t)$ by dividing the tumor cell population at time $t$ by the total number of grid positions $10,000$.

(3) Cell Stickiness $(S_i)$ - Cells can secrete substances that allow cells to stick with neighboring cells. The stickiness of cell $i$ is represented by $S_i \in [0, 1]$. In cell proliferation, healthy daughter cells take the same stickiness as their mother cell. Otherwise, the stickiness of the daughter cell differs from the mother cell, depending on the change in stickiness. The movement of a cell depends on the stickiness of that cell. We assume that healthy cells have a high stickiness value ^[37].

(4) Change in Stickiness $(\mu_s)$ - When a cell mutates, the stickiness of that cell changes according to the stickiness of its mother cell. Assume that the change in stickiness is $\mu_s$ and the stickiness of the mother cell is $S$, then the stickiness of the daughter cell can take a uniformly distributed value between $S-\mu_s$ and $S+\mu_s$.

(5) Jump Radius $(J_i)$ - The jump radius of cell $i$ is represented by $J_i$. It measures cell movement and represents the number of steps a cell can move at a time. We assume that the movement direction of a cell should be North, South, East or West. The jump radius of a tumor cell changes according to the jump radius of its mother cell. We assume that the jump radius of the daughter cell can increase by one, decrease by one or take the same jump radius as the mother cell, with a minimum value of $1$.

(6) Change in Jump Radius $(\mu_j)$ - When a cell mutates, the jump radius of that cell changes according to the jump radius of its mother cell. Assume that the change in jump radius is $\mu_j$ and the jump radius of the mother cell is $J$, then the jump radius of the daughter cell should take a positive integer value chosen uniformly in the range $[J-\mu_j, J+\mu_j]$. In this model, we assume that $\mu_j = 1$.

(7) Maximum Divisions of a Healthy Cell $(K)$ - Healthy cells can divide only a finite number of times. Due to genetic changes, tumor cells have the potential to reproduce unlimited times ^[92]. However, they have a slight chance of dying after completing the maximum number of cell divisions.

(8) Number of Divisions of a Cell $(D_i (t))$ - The total number of divisions of cell $i$ until time $t$ is represented by $D_i (t)$. For a healthy cell, $D_i (t)\leq K$. Since tumor cells can proliferate uncontrollably, divisions of a tumor cell can take any positive integer value. However, tumor cells have a slight chance of dying when they complete $K$ cell divisions.

(9) Cell Age $(a_i(t))$ - The age of cell $i$ at time $t$ is represented by $a_i(t)$. At each time step, the age of each cell increases by one. When a cell reaches its reproduction age, it has a chance to proliferate. When a cell proliferates, the age of the mother cell and the daughter cell initialise to zero, and the number of divisions of that cell ($D_i (t)$) increases by one.

(10) Division Age $(L)$ - A cell can divide when the age of that cell is greater than the division age. The duration of the cell cycle differs from organism to organism and cell to cell. In this model, we assume that the division age of a cell is $20$.

(11) Competition Rate $(C_i)$ - The competition rate of cell $i$ is represented by $C_i$. Each cell has a competitive rate that depends on the metabolism of that cell (the metabolism level of a healthy cell decreases with each cell division) ^[99,100]. If the maximum number of cell divisions of a cell is $K$ and the number of remaining cell divisions of that cell is $n$, the competitive rate of that cell is $n/K$. The range of values for the competition rate of healthy cells is $[0, 1]$. Since tumor cells have super-competitive phenotypes with their genetic changes, we assume that tumor cells have highly competitive rates. Hence, we choose a random competitive rate between $[0.5, 1]$. If a tumor cell does not have enough space to proliferate, that cell competes with immediate neighboring cells. If a cell in the neighborhood has a lower competition rate than the middle cell, that cell dies, and the daughter cell occupies that grid position.

(12) Probability of Dividing $(P_d)$ - If the age of a cell reaches the cell division age, the probability of dividing determines whether that cell divides or not at that time. Since healthy cells can control their cell growth, they only proliferate when they get signals to execute. Tumor cells have a high probability of cell division because they sustain proliferative signalling and evade growth suppressors.

(13) Probability of Death $(P_a)$ - When a healthy cell completes its maximum number of cell divisions, that cell has a higher chance to die. Since tumor cells resist their death, the probability of death of tumor cells is small. When a cell dies, the grid position of that cell becomes empty.

(14) Obstacle Density $(R)$ - The ECM provides a physical barrier for cell movements through proteins. In this study, we include obstacles in the domain to represent ECM proteins. If the next chosen grid position of a cell is occupied by an obstacle, the cell stays at its current position until it gets another chance to move. We assume that tumor cells cannot destroy obstacles. Hence obstacle density does not change with time.

(15) Probability of Cell Movement $(P_m(i, t))$ - The probability of cell movement depends on the number of cells or obstacles in the Moore neighborhood (eight surrounding grid positions) and the stickiness of the cell. The probability of movement of cell $i$ at time $t$ is represented by $P_m(i, t)$. If the stickiness of a cell is $S$ and the number of cells and obstacles in the neighborhood is $n$, the probability of movement of that cell is $(1-S)^n$.

(16) Fence Density $(f(t))$ - The boundary of the domain has a fence that represents ECM components of the TME. Generally, cells cannot move through the surrounding membrane ^[16]. Hence we assume that initially, the domain is fully fenced. Once tumor cells reach the boundary, they attach to the fence components and secrete proteases to break down the fence components. This process makes gaps in the boundary, and cells can move through these gaps.

(17) ECM Break Down Probability $(P_b)$ - When a tumor cell attaches to the ECM, some tumor cells secrete a high level of MMPs to break down the ECM component. Therefore some tumor cells that reach the boundary escape from the domain by making gaps in the boundary with a certain probability $P_b$.

(18) Number of Exit Cells $(E(t))$ - After reaching the boundary, tumor cells exit the domain through gaps in the boundary. The number of Exit Cells ($E(t)$) represents the total number of tumor cells that escaped the domain from time $0$ to $t$.

(19) Average Stickiness $(\bar{S}(t))$ - The output variable $\bar{S}(t)$ calculates the average stickiness of tumor cells at time $t$.

(20) Average Jump Radius $(\bar{J}(t))$ - The output variable $\bar{J}(t)$ calculates the average of jump radii of tumor cells at time $t$.

A description about the flow chart that represents in Figure 5

If the cell survives or $k\leq K$, we should check the age of that cell $(a_i)$. If $a_i\geq L$, the division age, generate a uniformly distributed random number $r_2\in[0, 1]$. That cell can divide if $r_2\leq P_d$. In the cell division process, the mother cell stays at the same position, and the daughter cell takes an unoccupied neighboring grid. The daughter cell placement can take a North, North-East, East, South-East, South, South-West, West or North-West grid, the first empty grid position in the clockwise direction starting from the North. If all neighboring grid positions are occupied, the mother cell competes with surrounding cells, and the weakest cell in the neighborhood dies, and that position is taken by the daughter cell. Figure 6 represents possible grid positions of the daughter cell after proliferation and how the position of the daughter cell can be chosen according to the cells or obstacles in the neighborhood.

Any alive cell can move to a new grid position depending on the number of cells and obstacles in the neighborhood, and cell stickiness. We generate a uniformly distributed random number $r_3\in[0, 1]$ and compare it with $P_m(i)$ of that cell that depends on the number of cells in the neighborhood and their stickiness. If $r_3 \leq P_m(i)$, the cell moves to a new grid position. A tumor cell can move North, East, South or West. If $r_3 > P_m(i)$, the cell stays at the same grid position. Moreover, cell movement distance depends on the jump radius of that cell. The randomly selected new cell position should be an empty grid position.