1.Introduction

According to the WHO ¹, cancer is the second leading cause of death worldwide, and it is estimated that the number of new cases will increase in the coming years. There is enough evidence ^2,3 regarding the complexity of cancer. We mean that cancer is seen as a self-organizing nonlinear dynamic system, far from thermodynamic equilibrium ^4,5,6, exhibits a fractal structure that allows it to evade the action of the immune system and host cells ^7-12. Shows growth dynamics in the metastatic stage, deterministic chaos type, which confers high robustness, poor long-term prognosis, adaptability, and learning capacity ^13-16.

As we have shown in previous works ¹⁶, cancer can be seen as a complex network made up of cells that have lost their specialization and growth control, and that emerges through what we can call “biological phase transition” ^17-19.

Cancer evolution and emergence are known to occur through three distinct stages (avascular, vascular, and metastatic) ²⁰; researchers often concentrate their efforts on answering specific questions on each of these stages. In the avascular phase, the tumor grows to a state known as a dormant state ²¹, with microscopic characteristics, this process, occurs through a “second-order” phase transition ¹⁹ through either logistic or Gompertz dynamic equations equivalently ¹⁹, which corresponds to the fact that such growth occurs silently.

The avascular phase, dormant state ²², is a stable stationary state, taking on a spheroid shape of a few millimeters in diameter, whose histological examination shows three different concentric annular layers ^23,24. Proliferating cells ⁷ are found in the thinner outer layer with a rough edge, fractal ⁷. In the contiguous layer, typically three times thicker than the proliferation layer, inactive cells, quiescent cells, exhibit little or no proliferation. The innermost nucleus consists of necrotic remains ²⁵. The tumor can remain in this latent state for months or even years, reaching a size between 1-2 mm in diameter and cannot be detected on a macroscopic scale ²⁶; therefore, its detection is usually limited.

This work aims to generalize the previously proposed model of avascular tumor growth ¹⁹ with the inclusion of spatiotemporal behavior. The manuscript is organized as follows: In Sec. 2, a two-dimensional (2D) cellular automata model is presented for the spatiotemporal study of avascular tumor growth. Section 3 describes the experimental part: the spatiotemporal avascular growth, including the dynamical behavior, the fractal dimension, and the entropy production rate of the tumor patterns. Finally, some concluding remarks are presented.

2.A Two-Dimensional (2D) Cellular Automata Model for the Spatiotemporal of Avascular Tumor Growth

Mathematical models ²⁷ represent a language for formalizing the knowledge on living systems obtained in theoretical biology ^28,29,30. Basic models of tumor growth ^31,32 make the description of the principal regularities possible and are a useful guideline for cancer therapy, drug development, and clinical decision-making.

These models can be classified into two general groups: deterministic based on ordinary or partial differential equations (ODEs, PDEs) ³³ and stochastic models ^34-37. On the other hand, to describe a spatial pattern formation has been used indistinctly, PDEs ³⁸, agent-based models (ABM) ^39,40, particularly, cellular automata (CA) ^41-43, hybrid cellular automaton ⁴⁴, and multi-scale models ⁴⁵.

CAs have been widely used given their ability to generate a broad spectrum of complex patterns from relatively simple rules which capture many behavioral characteristics and the complexity of real self-organizing systems far from thermodynamic equilibrium ³⁹.

The advantage of CA models is that they can directly incorporate the rules that define the processes of mitosis and apoptosis, as pointed out by Basanta et al. ⁴⁶ based on the presence of cancer cell hallmarks ⁴⁷.

The proposed CA cellular automaton is a reinterpretation of the model of Pourhasanzade et al.⁴⁸, a two-dimensional (2D) stochastic agent-based model for the spatiotemporal study of avascular tumor growth. For the development of the CA model and spatial patterns, the Python 3 programming language was used ⁴⁹.

The main differences with the model of Pourhasanzade et. al. ⁴⁸ are:

The proposed model is made up of a two-dimensional square lattice (2D, n x n cells); each cell in the lattice is defined by position coordinates (i,j), where 0<i,j≤n cells. These agents can take 4 possible modes: 0 - normal cells (N); 1 - proliferative cells (PT); 2 - non-proliferative cells (NT); 3 - necrotic cells (Ne). Besides the empty spaces in the lattice are shown as agents in mode 0.

The model includes a set of update rules at the cellular level of the biological system consisting of two sets of transition rules for agents. Thus, it considers proliferative cell mitosis (PT), competition between cell populations, and interaction between normal and tumor cells.

The transition rules for this model are probabilistic, and a Moore neighborhood [50] limits the interaction of each lattice site with its neighborhood. On the other hand, the nutrients are supposed to be evenly distributed in the lattice, where its deficiency can be considered due to the lack of free space within a certain distance from a cell. Therefore, a PT cell can only be divided with some probability into two daughter cells as long as there is a space in its vicinity, placing a daughter cell. In contrast, the other daughter cell will replace the position from which it originated.

Each mode 1 agent (PT cell) can proliferate with a probability as a function of time and space, such that

where P₀ is a base probability of proliferation, which in the original work by Pourhasanzade et al. ⁴⁸ it is assigned an arbitrary value, in our case, according to the temporal model previously developed by us ¹⁹ we associate by ansatz to the quotient of the mitosis (ψ) /apoptosis (η) constants, P0(ψ/η), whose quotient characterizes a given tumor ⁵¹. Consequently, the probability P of proliferative cells (PT) is associated with the rates of mitosis (V _m ) and apoptosis (V _a ), respectively.

The n₀ reflects the number of healthy cells (N) in each site’s vicinity, associated with the morphological characteristics of the host tissue. The r is the position of each cell with respect to the center of the tumor and R _max is a constant value that is prefixed and denotes the maximum diameter between 1-2 mm that the tumor can reach ⁵² since it is known that they grow in the avascular stage until they reach a latent state. Appendix A shows a flowchart showing the sets of rules for the agent-based model proposed.

In this way, it is ensured that the CA model results are adjusted through either logistic or Gompertz dynamic equations equivalently ¹⁹ and that it follows a dynamic where for radius greater than R _max the probability of proliferation is zero, thus showing the effect of the mechanic pressure exerted by the host on tumor growth ²⁴.

In each time iteration, if the agent mode is 1, it is checked if the agent has opportunities to proliferate. For this, there must be at least one agent with mode 0 in its vicinity. If so, the corresponding PT cell will choose one of those possible positions at random and divide with a probability P. Of the two new cells created, one takes the position of the cell in mode 0, and the other remains in the position of the PT cell.

In case the PT cell cannot proliferate, either because it cannot find a healthy neighbor in its vicinity or because it cannot proliferate with a probability P, it will remain as a PT cell for a certain time. Therefore, each PT cell has an age counter that is incremented at each time step and is reset if the cell undergoes mitosis. After reaching a certain time (age threshold), the PT cell can change to the NT cell, changing the mode from 1 to 2 when it is at a greater distance than δP, the one associated with a certain distance from the edge of the tumor.

NT cells that are in mode 2 can change to Ne cells (necrotic cells, mode 3) if they are at a greater distance than δP+δn from the edge of the tumor due to the absence of nutrients. These Ne cells are found in the center of the tumor as a mass of dead cells and do not change their mode under any circumstances.

The radius of the necrotic region (R _n ) and the thickness of the region of PT (δP) and NT cells (δn) are obtained by the following equations ⁴⁰:

where α and b are constant parameters that reflect the need for nutrients in tumor growth, and R_t is the average radius of the tumor. Table I shows the values of the parameters used in the CA model. In all cases, the parameters used are dimensionless.

3.Experimental Part: the Spatiotemporal Avascular Tumor Growth

Figure 1 shows the patterns that are obtained for different values P₀ after 200 iteration steps, when saturation is achieved, that is, the dormant state. In purple, it represents healthy tissue (N), blue is proliferative cells (PT), green is non-proliferative cells (NT), and yellow is necrotic (Ne).

Figure 1 Patterns formed during avascular growth obtained for different values of p ₀ (0.1; 0.3; 0.5; 0.7) after 200 iteration steps have reached the dormant state: healthy tissue (N, purple), proliferative cells (PT, blue), non-proliferative cells (NT, green), necrotic (Ne, yellow). 

For each fractal pattern (see Fig. 1), the fractal dimension d _f was determined using the box-counting method ^53,54. Each image was processed with the ImageJ 1.40 g software by Wayne Rasband, National Institute of Health, USA ().

As shown in previous work ³⁶, the fractal dimension d _f can be given as a function of the quotient between mitosis (V _m ) and apoptosis (V _a ) rates ^36,55, different values of d _f representing the “degree of malignancy” ⁵⁶ which quantifies the capacity of the tumor to invade and infiltrate the healthy tissue, as

From the Eq. (3) the quotient of mitosis/apoptosis rates (V _m /V _a ) was determined (see Table II), which increases as the value of growths P₀, which physically implies an increase in the degree of malignancy of the patterns formed. Notably, the fractal dimension values for the different patterns are in the range of values reported by Brú ⁷ for different tumor cell lines (see Fig. 1, Table II).

Figure 2 shows how the cell population varies with time, for different values P₀ after 200 iteration steps, so that saturation is achieved, that is, the dormant state (see Appendix B).

Figure 2 Temporal variation of the cell population (PT) for different values of P ₀ (0.1; 0.3; 0.5; 0.7) after 200 iteration steps. 

As observed, a sigma-shaped curve the avascular growth of the tumor, follows a growth dynamic, through either logistic or Gompertz dynamic equations equivalently ^11,19. As we pointed out in Sec. 2, the probability p of proliferative cells (PT) is associated with the rates of mitosis (V _m ) and apoptosis (V _a ), respectively. This shows how, for different types of tumors, characterized by their P_O value, the tumor exhibit varying degrees of aggressiveness ¹⁹. We previously reported ³⁶ the same behavior that in this work is obtained for different types of tumor cell lines.

Previously we have shown ⁵¹ that the entropy production rate Ṡi was determined for avascular tumor growth as a physical function of cancer robustness, such as

The Eq. (4) shows two major properties associated with the avascular tumor growth: The first is its growth rate, which is associated with its invasive capacity, mitosis (V _m ), and apoptosis (V _a ) rates, related with the degree of aggressiveness. The second is its complexity, a morphology characteristic, such as the fractal dimension of the tumor interface, associated with malignancy, which quantifies the tumor capacity to invade and infiltrate the healthy tissue ⁵⁶. Figure 3 shows the entropy production rate Ṡi for different values of p0.

Figure 3 The entropy production rate S˙i for different values of p ₀. 

As can be seen, the entropy production rate Ṡi shows how the robustness of the tumor increases as its invasive capacity increases, degree of aggressiveness, and in turn, its malignancy.

4.Conclusions and Remarks

In summary, in this paper, we have found that the proposed two-dimensional (2D) cellular automata model generalized the spatiotemporal behavior of avascular tumor growth, and that allows a better understanding of the morphogenesis of the tumor pattern formation.

Using thermodynamics formalism of irreversible processes and complex systems theory, we propose and quantify markers able to establish, in a quantitative way, the degree of aggressiveness and the malignancy of tumor patterns, such as:

The current theoretical framework will hopefully provide a better understanding of cancer growth and contribute to cancer treatment improvements.