Agent-Based Model for Studying the Effects of Solid Stress and Nutrient Supply on Tumor Growth

Abstract

An off-lattice agent-based model of tumor growth is presented, which describes a tumor as a network of proliferating cells, whose dynamics depend on the stress generated by intercellular bonds. A numerical method is introduced that ensures the smooth dynamics of the cell network and allows for relative numerical cheapness while reproducing the effects typical of more complex approaches such as the elongation of cells toward low-pressure regions and their tendency to maximize the contact area. Simulations of free tumor growth, restricted only by the stress generated within the tumor, demonstrate the influence of the tissue hydraulic conductivity and strength of cell–cell interactions on tumor shape and growth rate. Simulations of compact tumor growth within normal tissue show that strong interaction between tumor cells is a major factor limiting tumor growth. Moreover, the effects of normal tissue size and strength of normal cell interactions on tumor growth are ambiguous and depend on the value of tissue hydraulic conductivity. Simulations of tumor growth in normal tissue with the account of nutrients yield different growth regimes, including growth without saturation for at least several years with the formation of large necrotic cores in cases of low tissue hydraulic conductivity and sufficiently high nutrient supply, which qualitatively correlates with known clinical data.

1. Introduction

From the point of view of the basic unit of life—the cell—the evolutionary struggle results in radically different optimal strategies in unicellular and multicellular organisms. In the latter, cells have to cooperate with each other according to a rigorous program embedded in their genome and are controlled by intercellular signaling pathways [1]. This cooperation even includes a cell’s readiness to sacrifice itself for the good of the whole organism [2]. However, if the cell represents a whole organism by itself, its optimal strategy is the opposite and consists of staying alive by any means necessary and proliferating as quickly as possible in order to distribute its genes. Although cancers can emerge in multicellular organisms, the malignant cells prefer to adhere to unicellular organisms, which requires hacking the cooperation program and relying more on primitive selfish behavior patterns [3]. Notably, all the established distinguishing features of cancer are focused on increasing its proliferation rate either directly or by adjusting the conditions of its microenvironment [4,5].

A cancer cell can thus proliferate almost without limits under favorable conditions. However, for a tumor growing in living tissue, such conditions cannot be achieved by all cancer cells, except in the very early stages of tumor growth only. As cancer cells require material, energy, and free space for their proliferation, nutrient deficiency and elevated stress are the main factors than limit tumor growth. Although cancer cells have the mechanisms to counteract these restrictions, their effect is inevitably limited and generally supports only a fraction of malignant cells. For example, cancer cells generally have low sensitivity to the mechanical stress generated due to their proximity with other cells, which, in a normal case, should prevent the tissue from excessive proliferation [6]. However, local stress is not alleviated fast enough, e.g., through tissue displacement, a proliferating cancer cell ultimately reaches the point where its act of division is impossible since its cytoskeleton cannot withstand unlimited load during mitosis [7]. Therefore, the increase in solid stress as a result of the proliferation of tumor cells further inhibits their proliferation rate [8].

In order to overcome nutrient deficiency, cancer cells stimulate angiogenesis, i.e., the formation of new capillaries that supply them with nutrients [9]. However, as a rule, new capillaries are concentrated only near the tumor rim, whereas in deeper parts of the tumor, the capillary network is usually poor due to the collapse and rupture of microvessels caused by elevated stress [10]. The invasion of tumor cells into the surrounding tissue allows them to move away from the nutrient-deprived zones but obviously helps only a fraction of cells that have acquired this ability.

Solid stress and nutrient availability in tumors, therefore, are interesting from a therapeutic perspective. Antiangiogenic therapy has been used in clinics since 2004 [11] and various mechanotherapeutic approaches are gaining popularity in preclinical oncological studies [12,13]. The effects of these types of treatments, however, can be difficult to predict and vary widely. A prominent example is the clinically detected stimulation of tumor progression toward increasingly invasive phenotypes through antiangiogenic therapy [14]. Moreover, the efficiency of antiangiogenic therapy varies widely for the different types of tumors [11].

A relatively cheap and fast way to obtain insights into the reasons for the behaviors of different tumors and decipher the mechanisms of their reaction to various external factors is the investigation of mathematical models. In these models, a tumor and its microenvironment are considered a single complex system. This approach allows for not only an explanation of the observed effects [15,16] but also ideas for the direction of further experimental and clinical research [17,18].

There are two ways to mathematically model tumor growth when accounting for solid stress [19]. The first is continuous modeling, which involves partial differential equations to govern the model. In relatively simple cases, related works consider tissue a liquid-like medium [20,21]. More complex methods are adapted from solid mechanics [22,23]. However, such methods are associated with significant theoretical difficulties, as living tissues, especially those susceptible to significant deformations, are very different from non-living solids. Tumor cells divide and die, and the network of interconnected cells, as well as the extracellular matrix, undergo constant remodeling in response to deformations, which, in turn, affects the stress distribution [24]. Notably, the consideration of nutrient dynamics does not present significant difficulties and the use of reaction–diffusion equations for this purpose represents a natural, physiologically-based method.

The second way of mathematical modeling is discrete modeling, wherein the dynamics of each cell are considered separately, in particular, its position and speed are monitored. With this formulation, the behavior of the network of interconnected cells can be simulated quite naturally. However, these methods are typically significantly more computationally expensive than continuous modeling methods, in which the cell density and other observable parameters are averaged in space. However, discrete approaches have recently become very popular due to the rapid development of computer technologies. Notably, when nutrients are considered explicitly in such models, their dynamics are almost always still modeled using continuous reaction–diffusion equations.

Various discrete methods of varying complexities exist that have been applied to solutions of a broad range of tasks related to oncology (see, e.g., [25,26] for review). The simplest and most popular approach is a cellular automaton, where a regular, generally non-deformable, spatial lattice is introduced, of which each element can be occupied by a cell [27]. In the more complex, cellular Potts model approach, each cell occupies several elements of a regular lattice [28], which allows for the consideration of explicitly varying cell shapes but inevitably requires much more computational power.

Since the constraint of fixed-space lattice significantly limits the consideration of mechanical effects, the most popular approach for the modeling of tumor growth when accounting for solid stress is the use of off-lattice agent-based models in which continuous space is considered (see the reviews in [29,30]). The most frequent basic assumptions of this type of modeling are as follows. The cells are represented by elastic spheres or circles, which exert a certain force on nearby cells. Two cells push each other apart if they are too close and attract each other if they are far enough, with the interaction force vanishing at a certain distance. The forces applied to each cell cause the cell to move in the direction of its sum. The cell network structure undergoes remodeling due to cell movement, with new bonds forming between close cells and old bonds disappearing between distant cells. Such an approach can reproduce basic experimental phenomena such as the transition from exponential to sub-exponential growth at sufficiently large tumor sizes [31], as well as more specific phenomena such as the infiltration of inert microspheres into a tumor spheroid via advection from its surface [32]. Thorough calibration of such models enables obtaining quantitative correspondence between the model behavior and experimental data dependent on parameter variations (for cell numbers up to ${10}^6$) [33]. In particular, such calibration reproduces the quantitative response of cell population dynamics to the variation of externally applied mechanical stress [34]. Furthermore, such models allow for the proposal of new hypotheses at a qualitative level. For example, the work in [35] suggests that an increase in tissue stiffness can promote the aggressive morphology of tumor invasion despite slowing down the overall tumor proliferation rate, which is similar to the possible outcomes of antiangiogenic therapy mentioned above.

The key problems with such models are the descriptions of the processes of cell growth and cell division, as well as the description of the mechanism of influence of mechanical stress on these processes. One simple solution to these problems describes cell growth via the increases in its radius and division via its instantaneous transformation into two smaller round cells, which are displaced either in a direction where there is enough free space [36] or merely in a randomly chosen direction, implying that cells can overlap [32,33]. In the first case, local stress influences cell division only implicitly, whereas in the second case, it does not affect it. Obviously, such a description of cell division infers an abrupt change in the shape of a dividing cell and, therefore, an abrupt change in the interaction forces between it and the cells with which it is connected to. A way to overcome such abrupt changes is to model cell division via duplication of a cell, with two daughter cells immediately replacing it and almost completely overlapping at first [35]. Technically, due to the already introduced law of repulsion, this should enact a great repulsion force between these new cells, which should lead to a fast cell divergence that yields strong local deformations. These deformations could be smoothened by restricting the value of the repulsion force, but under significant local stress that counteracts the divergence of new cells, this would lead to the practical inability of cells to completely divide. A smooth process of cell growth can be ensured by regarding it as a separate process with specific laws. However, abandoning the typical spherical shape of cells and considering their more complex geometries significantly increases computational costs [37]. A computationally cheap solution would be to consider a cell as two overlapped spheres or circles that diverge at a moderate speed, corresponding to the typical tumor cell proliferation rate if the external stress is not too high [31]. However, this approach requires modifying the laws of connection between cells, as a direct account of the connections between each pair of spheres or circles is unnatural in this setting.

In this work, an off-lattice agent-based model of tumor growth is introduced, which includes a modification to ensure the smooth dynamics of the cell network while maintaining relative computational efficiency. We use this model to simulate tumor growth under three settings and explore its responses to variations in the biomechanical parameters. The first setting involves free tumor growth, which is only limited by the residual stress generated within the tumor due to the interconnection of tumor cells. In real tumors, this stress remains, even after the excision of tumors from tissue [38]. The second setting involves tumor growth within normal tissue, which adds the effect of external stress on tumor growth. The third setting considers the dynamics of a nutrient whose local level affects the behavior of tumor cells.

The computational code was implemented in C++ and can be found in the Supplementary Materials.

2. Free Tumor Growth

In this section, the model of a two-dimensional tumor growing freely, i.e., in the absence of normal tissue, is introduced and the results of its simulations are presented and discussed. As normal tissue is not yet considered, nutrient dynamics is also not taken into account and it is assumed that the behavior of tumor cells depends solely on the solid stress applied to them. Such a setting may correspond roughly to the growth of a tumor cell colony in a Petri dish with constantly replenished nutrients.

2.1. Model

2.1.1. Major Assumptions

The most basic model assumptions are quite general for this type of agent-based modeling and were described in Section 1. Briefly, cells are represented by circles, which push each other apart when too close and attract each other when not too distant from each other. Emerging forces cause cell movement and remodeling of the cell network structure. In the absence of cell division and death, the above description covers all the main features of cell dynamics on a qualitative level. However, for proliferating tumor cells, additional assumptions are introduced into the considered model.

In this work, a homogeneous tumor cell population is considered, where all tumor cells are identical and obey the same laws with the same set of corresponding parameters. Each of the circles that represent tumor cells never exists alone in this model—it is always paired with another circle and together they represent a growing and dividing tumor cell. This combination of two circles is referred to as a tumor cell pair, whereas each of the individual circles is still referred to as a tumor cell. Please note that this terminology is used for the convenience of the formulations and does not strictly correspond to the underlying biological foundation. In addition, please also note that the term “couple of cells” is used to refer to the unpaired cells, i.e., the ones that do not form a dividing tumor cell but have to be considered together. The paired tumor cells always overlap, i.e., the distance between their centers is lower than the sum of their radii. In the absence of external forces, the paired cells diverge from each other at a constant speed. When they stop overlapping, they become connected (and at first, repulse each other) and each of the cells becomes paired with a newly introduced one whose position is just a tiny bit shifted from its own. The forces applied to the paired cells due to the connections with other cells can slow down their divergence to a complete stop. Additionally, due to the external forces, paired cells move in space, with this movement being a combination of parallel transfer and rotation with respect to their center of mass. However, this method of considering tumor cells creates certain difficulties in defining bonds between cells and calculating the corresponding forces, which are explained below together with a solution.

In real tumors, solid stress is exerted on cells not only by other cells but also by the extracellular matrix. In this model, the extracellular matrix is not considered explicitly. However, implicitly, its density is assumed to influence the value of tissue hydraulic conductivity and thus the speed of cell displacement. Moreover, in Section 3, restricted tumor growth is considered, corresponding to the growth of benign tumors in capsules consisting of extracellular matrix elements.

2.1.2. Mathematical Formulation

The model simulations considered in this section represent a cyclic sequence of the following stages, with each cycle happening during a specific time period $\tau$. The values of the parameters introduced below are provided in Section 2.1.3.

The first stage is cell movement. Each pair of tumor cells participates in three separate types of movement: cell divergence, rotation, and parallel transfer. The divergence of cells mimics the increase in tumor cell volume before mitosis and, eventually, this process results in overall tumor growth. Each of them is influenced by the forces external to the considered pair of cells. The specific way of calculating these forces is explained below, but suffice it to say that they are directed along the lines that connect the centers of cells (recall that each of the circles is referred to as a cell) and the value of each force is a product of several coefficients, whose own product is positive, and the stress $\sigma$, whose value depends on the length of bond l. In considering free tumor growth, all the cell radii r are constant and equal to $r_0$. The following function of stress is used, which is adapted from the works that use continuous models [21,22] and is illustrated in Figure 1.

$$\begin{array}{c}\sigma \left(l\right)=k_T[c-c_0]{[c-c_{br}]}^2/\sqrt{c_{inf}-c},c={\{2r_0/[2r_0+l]\}}^3\\ c_0=c\left(l_0\right),c_{br}=c\left(l_{max}\right),c_{inf}=c(-l_0).\end{array}$$

(1)

This function assumes that at the normal bond length, $l=l_0$, cell interaction yields zero stress. If cells are pushed together, $l<l_0$ and a repulsive interaction appears. If cells are pulled apart, $l>l_0$ and an adhesive interaction at first strengthens and then weakens with the increase in the bond length, which corresponds to the sequential ruptures of individual intercellular contacts. For sufficiently separated cells at $l>l_{max}$, there is no interaction between them. Parameter c roughly corresponds to the volume fraction that cells would locally occupy in a three-dimensional space under all cell bonds equal to l in that location (neglecting the specific form of cells, i.e., assuming that under $l=0$, the whole volume is occupied by cells). The bond length can be negative; thus, the non-paired cells can overlap, which may be related to cell compressibility. However, as seen in Figure 1, this overlapping is not encouraged. Note that as cells approach each other, the repulsive response takes place before direct cell–cell contact ($l=0$), which reflects the fact that it is caused not only by the elastic deformation of cell membranes but also by the interaction of receptors and ligands attached to cell membranes [39].

Cell divergence. In the absence of external forces, paired cells move apart with a constant speed of $v_{div}=v$ along the line that connects their centers. When external forces are present (which was almost always true during the calculations), cell divergence slows down. It is assumed that cells cannot diverge under large enough pressure applied against the direction of cell displacement. That is, if an outer cell is bound with one of the considered paired cells in such a way that the centers of the three cells are on a single line and the outer bond length yields a pressure that exceeds ${\sigma }_{cr}$, then the considered paired cell cannot be displaced (its pair is, however, free to move apart if no external forces are applied to it). In general cases, the speed of divergence of cell X is

$$v_{div}^X=v·min(1,max(0,1-{\sigma }_D^X/{\sigma }_{cr})),$$

(2)

which is a non-negative value not exceeding v. Therefore, during each cycle, each cell X is moved using a numerical algorithm by the distance $d{R}_{div}^X=v_{div}^X\tau$. Here, ${\sigma }_D^X$ is the sum of the projections of all forces acting on cell X on the vector starting at its center and pointing at the center of its pair x, divided by the normalized area of the surface of interaction between two cells. The definition of this surface area in a two-dimensional case is a non-obvious task in itself so it is only accounted for implicitly. It is assumed that it depends only on the radii of the interacting cells, as discussed below, whereas its alterations due to the changes in the intercellular distance are not accounted for as they should bring only quantitative corrections to the stress function when they are multiplied. Since the stress function itself is of only a qualitative nature, such a refinement would be excessive. For two cells with radii $r_0$, however, the interaction surface area is considered to be normalized to unity. As cells diverge, the effective fractional area $f_S^X=f_S^x$ of each cell is monitored, which is defined so that $2f_S^X\pi r^2\in (\pi r^2,2\pi r^2)$ is the area occupied on the plane by the overlapped paired cells. It can be calculated from simple geometric reasoning that

$$\begin{array}{c}{\displaystyle f_S=1+\frac{sin\alpha -\alpha }{2\pi },\alpha =2arccos\frac{D}{2r},}\end{array}$$

(3)

where D is the distance between cell centers (superscripts denoting cells are omitted here and elsewhere for brevity).

Cell transfer. The external forces applied to a pair of cells X and x cause their parallel transfer, which is expressed in a traditional form reflecting the friction-dominated overdamped motion. This component of displacement is the same for both paired cells, which move along the direction of the sum of external forces applied to both of them, ${\mathbf{F}}^{\mathbf{X},\mathbf{x}}$, and during each cycle, they pass the distance

$$\Delta R_{tr}=K{F}^{X,x}\tau /\left[2f_S\pi r^2\right],$$

(4)

where the coefficient K is referred to as the tissue hydraulic conductivity, or tissue conductivity for simplicity. The denominator corresponds to the mass of a tumor cell pair, reflecting the assumption that greater cells are harder to move. The specific cell mass is not calculated explicitly but it is assumed that it should be proportional to the cell area. It can be speculated that the corresponding coefficient of proportionality is hidden in K since the tumor mass is only used in this model in similar fractions that involve tissue conductivity.

Cell rotation. During each time step, the paired cells rotate around their center of mass, which is located in the middle of the line segment connecting their centers, with the rotation angle being

$${\displaystyle \Delta \varphi =KT\tau /I,I=\frac{r^2}{4}\left\{2[D^2+2r^2]\left[\pi -arccos\frac{D}{2r}\right]+3D\sqrt{4r^2-D^2}\right\}}$$

(5)

Here, T represents the total torque about the cells’ center of mass, which is calculated directly, and I denotes the corresponding moment of inertia (note that the cell mass is implicitly accounted for here). It can be proven that the same coefficient K should be used in this formula by considering a limit case of a point body rotating around an infinitely distant axis, where the formulas for cell rotation and parallel transfer should yield the same results.

Redefinition of bonds. After cells have changed their positions, the distances between them change, which influences the length of each bond l and its existence or non-existence due to breaking. The most direct algorithmic way to track the bonds would involve checking each couple of non-paired cells during every cycle. However, a large number of cells would lead to a colossal number of operations during each time step (proportional to the square of the number of cells) and a vast amount of memory. In order to optimize the bond tracking, a less direct approach is used here, where operations with bonds are divided into the stages of bond redefinition (which includes their breaking) and new bond formation. The former deals with already existing bonds, the list of which can be found at the beginning of this stage.

As mentioned above, the consideration of paired tumor cells creates certain difficulties, which are considered now using the method illustrated in Figure 2, where dark gray circles represent paired tumor cells and light gray circles represent single normal cells (not present under free tumor growth but useful to consider). During the stage of redefining the bonds, the program goes through existing intercellular bonds and modifies the information about them. Each of the bonds is identified by two indices of interacting cells X and Y. Each bond is a line segment lying on the line that connects the centers of interacting cells and its ends are situated on the circumferences of the circles that represent these cells. The distance between the bond’s ends in a general case is modified during this stage, as $D^{X,Y}-r_X-r_Y$, where $D^{X,Y}$ is the actual distance between cell centers and $r_X$, $r_Y$ are their actual radii. If the bond does not go through cells that are paired with X and Y (see, e.g., the bond between tumor cell 6 and normal cell 8 in Figure 2), the length of the bond $l^{X,Y}$ is taken to be equal to the distance between the bond’s ends. These bonds are referred to here as “open bonds”. An important assumption of the model is that a tumor cell can interact with an outer cell with which its paired cell already interacts so tumor cell 7 and normal cell 8 are also connected. Note that if only the closest cell in a pair could interact with an outer cell, this could lead to abrupt reconnections during simulations. Such an approach could yield numerical difficulties when the whole cell network is considered, such as artificial oscillations of cells and rapid changes in the directions of cell rotation. The approach suggested here aims to avoid such abrupt changes and ensure the smooth dynamics of the cell network.

If only cells 6, 7, and 8 are considered, it is clear that the bonds with normal cell 8 should hold tumor cells in a position where the line connecting the centers of tumor cells 6 and 7 is perpendicular to the line connecting their center of mass with the center of cell 8. Such a feature corresponds to maximizing the contact area between cells, which can be explicitly implemented in cellular Potts models (see, e.g., [40]). In order to account for the increase in the cell contact area during the divergence of the tumor cells, the formula for the force generated by each bond accounts for the effective tumor cell fractional area $f_S$. For a general case, the force generated by a bond between cells X and Y is calculated as follows:

$$\begin{array}{c}{F}^{X,Y}=\sigma \left(l^{X,Y}\right)·k_l^{X,Y}·f_S^X·f_S^Y·min(r_X,r_Y)/r_0,\end{array}$$

(6)

where $\sigma \left(l^{X,Y}\right)$ is the stress that was defined in Equation (1) and $k_l^{X,Y}$ is the coefficient of the bond strength that is equal to 1 for open bonds. The last multiplier represents the implicit account of the area of the cell interaction surface. It accounts for the fact that smaller cells should have a smaller contact area and is equal to 1 when all cells have equal radii of $r_0$ (which is true during free tumor growth, as noted above). The effective fractional area of normal cells is always equal to 1. Therefore, when tumor cells 6 and 7 have just begun to diverge and thus their positions practically coincide in space, they exert a total force equal to $\sigma \left(l^{X,Y}\right)$ on cell 8, which is identical to the force that would be generated by a single normal cell in the same position under the same solid stress function. When the tumor cells are ready to separate, each of them acts with the same force on cell 8, provided that the bond lengths are kept the same.

If the bond between cells X and Y goes through at least one of the cells x and y that are paired with X and Y, this bond is considered a “closed bond” and its length is $l^{X,Y}=D^{X,Y}-r_X-r_Y-L_{Hx}^{X,Y}-L_{Hy}^{X,Y}$, where the last two terms represent the lengths of the parts of this bond that lie within x and y. They are referred to here as “hidden bond parts” and can be calculated through direct geometric reasoning. The hidden parts are marked with dashed lines in Figure 2, where the solid lines correspond to the “visible parts” whose lengths are equal to the defined bond lengths.

Consider normal cell 1 and tumor cells 2 and 3 whose centers are located on the same line. Note that the open bond between cells 1 and 2 and the visible part of the bond between cells 1 and 3 coincide. Imagine that tumor cells diverge while the lengths of bonds $l^{1,2}=l^{1,3}$ remain constant. In order to produce the constant total force with which tumor cells act on cell 1 in this case, the closed bond between cells 1 and 3 is assigned a coefficient $k_l^{1,3}=[1-f_S^3]/f_S^3$. Thus, the sum of two corresponding forces, calculated using Equation (6), is maintained equal to $\sigma \left(l^{1,2}\right)$ during tumor cell divergence. A similar thought experiment for two tumor cell pairs with all cell centers situated on a single line enables the deduction that in the case of closed bonds with two hidden parts (which connect the cells that are farthest from each other), the coefficient of the bond strength should represent a product of an analogical fraction. In a general case of random cell alignment, in order to ensure smooth dynamics, it is clear that $k_l$ should tend to zero when the lengths of the hidden bond parts are close to zero and the bond is going to become open. These features lead to the following general formula for the coefficients of the bond strengths:

$$\begin{array}{c}{\displaystyle k_l^{X,Y}=[1-\frac{2f_S^X-1}{f_S^X}·\frac{L_{Hx}^{X,Y}}{D^{X,Y}}]·[1-\frac{2f_S^Y-1}{f_S^Y}·\frac{L_{Hy}^{X,Y}}{D^{X,Y}}],}\end{array}$$

(7)

which is equal to 1, particularly for open bonds, as noted above. Note that in Figure 2, the bond between cells 3 and 4 is actually closed, because its hidden part, which goes through cell 2, is too small to be discernible. Therefore, its coefficient of strength $k_l$ is close to unity. Conversely, bonds with two hidden parts (such as between cells 4 and 7 and between cells 2 and 7) have the smallest values of $k_l$. Note that when defining the speed of cell divergence (see Equation (2)), closed bonds are assumed to affect the cells whose circumference the ends of their open parts are located on, e.g., the compressed closed bond between cells 2 and 7 slows down the divergence of their pairs, cells 3 and 6.

It is assumed that the bonds break when their length $l>l_{max}$. Thus, before breaking, the force generated by a bond vanishes, which contributes to smooth model dynamics. A broken bond is not considered during the next stage of bond redefinition.

Division of cells. If paired cells X and x have moved sufficiently far apart so that $D\ge 2r_X=2r_x$ and $f_S^X=f_S^x=1$, they are no longer considered a pair. They immediately become connected with each other through a new bond and a new cell pair emerges for each of them. The cells in each of the newly formed pairs are displaced in the opposite direction by a small distance $\Delta R$ from the position of the original cell. The direction of displacement is defined as perpendicular to the direction of the strongest force acting on the original cell. Thus, the cell pairs aim to avoid elongating in the direction of maximum pressure. This feature is maintained due to the ability of cells to rotate in response to external forces and the fact that the forces acting perpendicular to the direction of cell divergence do not hinder it. For each of the new cells, new bonds are assigned with all the cells that were connected to its actual pair, as well as with the cells of the other pair generated as a result of this division. The features of the new bonds are defined as described above.

New bond formation. Along with the formation of new bonds during cell division, new bonds should form when two unpaired cells become sufficiently close to each other, which is checked in the last stage. Since the direct checking of each cell couple during each cycle is a quite cumbersome task, this process is organized in the program as follows.

At first, the notion of the “potential bond” is introduced. Two cells are assumed to be potentially connected if the length of a bond that would connect them (which is defined according to the rules for closed bonds described above) is less than a predefined technical parameter, $l<l_{pot}$. The existence of potential bonds is checked for each couple of cells that are not paired, not connected, and do not already have a common potential bond. This action is rarely performed, certainly not during each time step, with a period of $T_{pot}>\tau$. More often, with a period of $T_{bond}\gtrsim \tau$, these potential bonds are inspected and each of them ceases to exist if its length is greater than another parameter $l>l_{pot}^{br}>l_{pot}$ and it turns into a (real) bond if $l<l_{max}$. Note that some other value smaller than $l_{max}$ could be used as a threshold for new bond formation but this value was chosen to provide smooth dynamics by avoiding the leaps in the force between two cells from zero to a notable finite value upon their connection and to keep the stress function $\sigma \left(l\right)$ simple (in particular, avoiding hysteresis effects).

If this stage takes place during the cycle, it is completed and after the increment in the time counter $t=t+\tau$, the program proceeds to the next iteration of the cycle.

2.1.3. Parameters

The set of parameter values used during the simulations of the free tumor growth is listed in

Table 1. Model parameters used during the simulations of free tumor growth.

Parameter	Description	Model Value(s)
Main:
$r_0$	cell radius (yet unchangeable)	5
$k_T$	coefficient of stress function for tumor cells	80/200/600
$l_0$	normal bond length (yielding zero stress)	0.75
$l_{max}$	maximum bond length	5
${\sigma }_{cr}$	critical stress at which cell divergence stops	15
v	maximum speed of divergence of a paired tumor cell	≈0.22
K	tissue conductivity	10/0.32/0.01
Technical:
$\tau$	time step	$(0.0004-0.004)/K$
$\Delta R$	initial displacement of paired tumor cells	5 × 10${}^{-3}$
$T_{pot}$	period of checking for potential bond formation	$0.28/K$
$T_{bond}$	period of checking for bond formation	$0.004/K$
$l_{pot}$	length at which potential bond is formed	25
$l_{pot}^{br}$	length at which potential bond is discarded	35

. Their dimensionless values were obtained using the following normalization values: 1 h for time, 1 $\mathsf{\mu }$m for length, and 1 kPa for stress. Note that for convenience, the values of the forces are not converted to some of their standard units in this work (as well as the values of masses, which have been discussed above). Direct estimations using the standard parameter values, however, suggest that in the three-dimensional case, the adhesion forces between two cells could realistically reach values of up to $0.5$ $\mathsf{\mu }$N (see, e.g., [41]).

The value for the tumor cell radius was chosen to provide a typical tumor cell size [42]. As stated in Section 2.1, the parameters of the stress function were based on the works in [21,22] that deal with continuous models. In accordance with the reasoning provided in our previous work [43], which also considered a continuous model, the value of the critical stress for cell proliferation was roughly assessed by the extrapolation of the data in [22]. The stress function coefficient $k_T$ varied depending on the functions, as depicted in Figure 1. Its greatest value led to the maximum absolute value of adhesion stress close to ${\sigma }_{cr}$, whereas at its minimum value, direct cell contact produced repulsion stress less than ${\sigma }_{cr}$, which stimulated cell overlapping during active proliferation. The value for the tumor cell pair divergence speed v was obtained as $r_0/[ln2/B]$, where $B=0.03$ was a reasonable value for the tumor cell proliferation rate, which was also used in our previous works (see [43,44]). A single tumor cell division unconstrained by external forces happened in about 23 h. The values of tissue conductivity also varied during the simulations, based on the assumption that since ${\sigma }_{cr}$ was the critical pressure that stopped tumor cell divergence, the following inequality must hold: $K{\sigma }_{cr}/\left[\pi r_0^2\right]\gtrsim v$. Tissue conductivity is known to vary significantly among tumors of various cell lines (see, e.g., [45]) and in this study, it varied across quite a broad range.

The initial displacement of paired tumor cells was chosen to be small compared to the cell radius—in the absence of external forces, tumor cells would pass such a distance in about 80 s. The time step was adjusted dynamically within the range depending on the value of K, avoiding significant alterations of the bond lengths, which could lead to numerical errors. The time step and values of the other technical parameters were adjusted empirically so as not to distort the solution or minimize the computational time spent on it.

2.2. Results

In this section, the results of the simulations of tumor growth in the absence of normal tissue are demonstrated. Tumor growth was restricted only by the residual stress that accumulated, as tumor cell numbers increased, along with the number of intercellular bonds. The external forces applied to the tumor were absent and no nutrients that would influence the behavior of the tumor cell were considered.

The initial conditions were a single pair of tumor cells, which were separated by a small distance of $2\Delta R$. Figure 3 shows five snapshots of the typical beginning of tumor growth, with this simulation performed under $k_T=200$ and $K=0.32$. At first, a single pair of tumor cells did not experience external pressure and divided in the minimum possible time, ≈23 h. After the division, two new tumor cell pairs formed immediately, oriented perpendicularly to its direction. Four intercellular bonds emerged. The crossed bonds were stretched and had a non-zero projection in the new direction of division but affected it only slightly. Compared to the first unconstrained division, it took ≈5% more time for each of the new cell pairs to divide. The two other bonds were compressed and compensated for the forces that tend to bring cells closer together but have no influence on the cell division rate.

During several first divisions, cells diverge in the direction perpendicular to the previous direction of their divergence, which should lead to a square arrangement of tumor cells. However, after several divisions, this type of arrangement breaks due to minor instabilities, and a more stable hexagonal (“honeycomb”) arrangement of cells manifests itself. At 97 h of tumor growth, this arrangement is rather accurate, with the lengths of the bonds differing from their normal stress-free length $l_0$ by no more than $17\%$. However, as proliferation continues, the cells located in the interior regions start to experience significant pressure, which slows down their division. As cells are interconnected and adhere to each other, such pressure cannot be readily alleviated by cell displacement. In the sixth generation of tumor cells, the process of division takes up to 34 h. Upon further proliferation, the division rate of cells located in the center of the tumor further slows down, whereas cells at its periphery divide at a rate close to maximum.

Figure 4 shows the dynamics of the number of cell pairs under free tumor growth with the variation of two parameters. The first parameter is the coefficient of the stress function for tumor cells $k_T$. Under an increase in its value, a couple of connected cells tend to both repel and attract each other with greater force when the length of their bond l deviates from its normal value $l_0$ (see Figure 1). The absolute values of maximum adhesive stress that a bond can generate under the three used values of $k_T$ constitute ≈13%, 33%, and 99% of the stress value, which stops tumor cell divergence when applied against its direction, ${\sigma }_{cr}$. Under the lowest value of $k_T$, the two cells that touch each other generate stress approximately equal to ${\sigma }_{cr}$, which stimulates a slight overlapping of non-paired cells during simulations. Under the highest value of $k_T$, the value of ${\sigma }_{cr}$ is reached under contraction of the bond down to ≈71% of its length. The second parameter is the tissue conductivity K, which correlates with the speed of cell displacement in response to the applied forces. Thus, a higher value of this parameter corresponds to a greater pressure alleviation within the tissue.

The highest rate of tumor growth was manifested under the highest of the used values of $K=10$ and the lowest value of $k_T=80$, when tumor cells were most readily displaced and could adhere only slightly to each other, yielding quite low resistance to displacement. As the corresponding graph shows, the number of tumor cells changed rather abruptly in similar periods, indicating that cells were sufficiently synchronized during the whole simulation. The best fit by an exponential curve suggests that in this case, the average cell proliferation rate fell by only about 8% compared to the unrestricted case, whereas by increasing $k_T$ to 200 and 600, it fell by ≈11% and ≈17% and cell synchronization was wiped out. At $K=0.32$, the growth curves became notably sub-exponential as the proliferation of interior cells was halted almost completely after a couple of weeks of tumor growth. The reduction in the overall tumor proliferation rate was more significant at $K=0.01$. In this case, after the initial 12 days, the best fit exponential curve for any $k_T$ corresponded to only about half the maximum tumor cell proliferation rate and quite a large number of cells effectively stopped proliferating.

Figure 5 shows the arrangements of tumor cells for the discussed simulations when the number of tumor cell pairs had just reached 3500. The colors of cells correspond to their speed of divergence $v_{div}$, with light gray cells dividing slower and black cells dividing at a speed close to maximum. This figure confirms that the overall tumor proliferation rate largely depended on tissue conductivity. Independent of the value of $k_T$, at $K=10$, actively dividing tumor cells were met throughout the tumor volume. At $K=0.32$, only the tumor rim actively proliferated and tumor cells dividing at close to the maximum speed were rarely met. At $K=0.01$, cell proliferation was significantly inhibited, even at the tumor rim since even the cells aligned along the tumor surface can experience significant stress in the direction of their divergence. Some of them escape it by passively turning in response to external forces and changing the direction of division but this process is rather slow due to the low conductivity. Strong cell–cell interactions that occurred at $k_T=600$ favor the clustering of cells yielding notable outgrowths whose widths decreased with the decrease in the tissue conductivity that suppressed cell proliferation within them. Lower values of the stress function coefficient, along with the moderate sizes of proliferating tumor rims, achieved at $K=0.32$ conversely led to more circular tumor shapes.

For each tumor, a graph is attached, which provides the distributions of the bond lengths for inner, middle, and outer cells, their positions defined with respect to the initial position of the center of the first cell pair. All these distributions are unimodal, with the peak values achieved at lengths smaller than $l_0$. Confirming the discussion in Section 2.1.3, the lowest value of $k_T=80$ stimulated cell overlapping under a significantly low tissue conductivity, even for the cells at the tumor rim. The increase in the value of $k_T$ and, therefore, stronger cell–cell interaction and a decrease in the value of K both led to the smoothing of differences between the distributions of bond lengths in different layers. With an increase in $k_T$, the main modes shifted close to $l_0$, as the cell–cell interactions became stiffer. With a decrease in K, more and more outer-located cells became arranged more compactly, leading to a fraction of the dividing cells eventually decreasing, even at the tumor border.

3. Encapsulated Tumor Growth in Normal Tissue

In this section, the model is augmented with normal tissue, the main element of which is a network of interconnected normal cells. In order to roughly simulate the elastic response of normal tissue to tumor growth, it is situated within the round area whose boundary is referred to here as a capsule. The capsule can be enlarged when sufficient pressure is applied to it from within, whereas by itself, it tends to shrink. This method also prevents excessive disruption of the network of interconnected normal cells in response to tumor growth. In the simulations in this section, the tumor is situated within its own round capsule, which separates it from normal cells. This capsule corresponds to real capsules made from extracellular matrix elements that can surround tumors, inhibiting tumor cell invasion in normal tissue and slowing down the tumor proliferation rate [46]. The newly introduced elements—capsules and normal cells—exert external stress on the growing tumor in addition to the residual stress discussed above. Furthermore, nutrients are introduced into the model, with the behavior of tumor cells dependent on their local level. For simplicity, it is assumed that the positions of capillaries should correlate with the positions of normal cells, and, therefore, the latter formally act as the sources of nutrients in the model.

3.1. Model Extension

3.1.1. Normal Cells and Capsules

The description of normal cells and their dynamics is analogical to that of tumor cells, which was described in Section 2.1.2, with the significant simplification that normal cells do not divide and each of them exists alone. The radius of each normal cell is always equal to $r_0$. As normal cells are considered round, their rotation is neglected and the only component of their motion is the parallel transfer governed by the law analogical to Equation (4) but with $\pi r_0^2$ in the denominator. They adhere to and repel each other by the law analogical to Equation (1), with, in a general case, a different stress function coefficient $k_N$. All the above-described logic related to the intercellular bonds and forces also applies to the normal cells, with their effective fractional area $f_S$ always equal to one.

The initial conditions in the following simulations always include a population of normal cells arranged in a hexagonal arrangement with equal stress-free bond lengths $l_0$ (as Figure 3 shows, tumor cells also tend to such an arrangement naturally). Normal tissue is encircled by a normal capsule with an initial radius $r_N=r_N^0$ and zero thickness, whose center should coincide with that of one of the normal cells and is considered the framework origin. Only the cells that fit in a normal capsule area are considered from the beginning of the simulation. The central cell is, however, replaced with a pair of tumor cells, displaced by a small distance of $2\Delta R$, and encircled by a tumor capsule, concentric with the normal one and with an initial radius of $r_T=r_0+l_0/2$ and a zero thickness. The first tumor cell division takes place immediately from the beginning of the simulation, yielding four overlapping tumor cells. This was performed in order to overcome the artificially caused inability of tumor cells to undergo the first division; otherwise, in the beginning, two dividing tumor cells have to counteract the pressure applied on the external surface of the tumor capsule via its bonds with six normal cells. In order to avoid the associated numerical difficulties, the solid stress function $\sigma \left(l\right)$ is modified here using the above restriction with a value of 10 $·{\sigma }_{cr}$, and is equal to it when $l\le -l_0$.

In the simulations, the main cycle, as described in Section 2.1.2, was preceded in each time step by two new stages that defined the dynamics of tumors and normal capsules. The centers of the capsules were held at the framework origin and the capsules maintained a circular shape. Therefore, only their radii $r_N$ and $r_T$ changed in response to the pressure exerted by cells. The bonds between cells and capsules were defined analogically to the intercellular bonds with the difference that cells can connect to capsules from within. In this case, a bond is well aligned along the line connecting the centers of a cell and a capsule. The length of the open bond, i.e., the one that does not go through the pair of a considered cell, is $l=r_C-r-R$, where $r_C$ is the radius of the corresponding capsule, $r_N$ or $r_T$, r is the cell radius and R is the position of its center. Here, closed bonds are defined analogically to the intercellular closed bonds and the forces between cells and capsules are calculated according to Equations (1) and (6), with, in a general case, a different stress function coefficient $k_C$. The parameter of the effective fractional area related to the capsules is always equal to 1. Under such conditions, tumor cells do not leave the region enclosed by the tumor capsule, therefore, compact tumor growth is simulated. Normal and tumor cells cannot connect through the tumor capsule.

During each cycle, the radii of the capsules $r_T$ and $r_N$ change according to the following values:

$$\begin{array}{c}{\displaystyle \Delta r_T=K_C^T\frac{F_T}{2\pi r_T}\tau ,\Delta r_N=K_C^N[\frac{F_N}{2\pi r_N}-s_0^N]\tau ,}\end{array}$$

(8)

where $F_T$ and $F_N$ are the sums of the forces acting on the corresponding capsules.

3.1.2. Nutrients and Modes of Cell Behavior

When nutrients are considered, they are represented by one generic nutrient for simplicity. It is supplied from the normal tissue (as stated above, normal cells effectively play the role of its sources), it diffuses within the tissue, and it is consumed by tumor and normal cells. Its local level influences the behavior of tumor cells. Calculations of the actual distribution of nutrient $n(x,y,t)$ represent an additional stage of the numerical algorithm performed at the beginning of every time step. Rather than using direct but numerically expensive finite-element methods that would be typically used in this situation, the following simple approach was used. Its use was justified by the fact that possible quantitative differences between methods are acceptable here as only the qualitative behavior of the tumor and normal tissue are of interest in this study.

The nutrient profile was calculated as if distributed on a network of interconnected nodes, each of which was situated at a normal cell center or a tumor cell pair center of mass. During each cycle, the following increments were added consecutively for nutrient levels at each of the nodes at the normal cell centers:

$${\displaystyle \Delta n_N=\{P[1-n]-Q_N\frac{n}{n+n^{*}}\}·{\tau }_n}$$

(9)

and at the tumor cell pair centers of mass:

$${\displaystyle \Delta n_T=-\left\{[{\gamma }_{pr}{Q}_T+Q_N]\frac{n}{n+n^{*}}\right\}·\frac{2f_S{r}^2}{r_0^2}·{\tau }_n,}$$

(10)

These terms represent the diffusive inflow from blood with coefficient P, with the nutrient level in blood normalized to unity, and consumption by tumor and normal cells, the latter obeying classical Michaelis–Menten kinetics. The nutrient consumption rate by a normal cell $Q_N$ was constant and tumor cells additionally needed a significant amount of nutrients for their proliferation; therefore, their nutrient consumption rate can increase up to $Q_T+Q_N$. The nature of coefficient ${\gamma }_{pr}\in [0,1]$ is explained below. Furthermore, the diffusion of nutrients was approximated at each node A by the addition of the following increments to the nutrient levels, where A corresponds to the index of a normal cell or one of the indices (the smaller one) of a tumor cell in each of the tumor cell pairs:

$${\displaystyle \Delta n_{dif}^A=\sum _{bonds+p.b.}\Delta n_{dif}^{AB},\Delta n_{dif}^{AB}=D_n·[n^B-n^A]·\frac{{\tau }_n}{{\left[d^{AB}\right]}^2},}$$

(11)

where the sum is taken across all nodes that correspond to cells that have bonds or potential bonds with cell A, $D_n$ is the nutrient diffusion coefficient, and $d^{AB}$ is the distance between nodes. The time step ${\tau }_n$ used in all three equations was adjusted automatically and this step was repeated N times, where ${\tau }_n=\tau /N$, in order to exclude numerical errors by keeping nutrient levels within the range $[0,1]$.

Under the consideration of tumor growth in tissue, the local nutrient level affected the behavior of tumor cells in a threshold manner so that cells could not only diverge but also shrink, die, and grow. Parallel transfer and rotation, governed by Equations (4) and (5), happened under all circumstances, but paired cells diverged only if the local nutrient level was higher than $n_q$ and cell radii were equal to $r_0$. It is assumed that for divergence, they need to consume nutrients at a rate proportional to their speed, therefore, for diverging cells X and x:

$${\displaystyle {\gamma }_{pr}={\gamma }_{pr}^{div}=\frac{v_{div}^X+v_{div}^x}{2v}.}$$

(12)

If $n>n_q$ but $r<r_0$, cells do not diverge but grow, i.e., their radii increase, with the center of mass position remaining fixed, by the following amount:

$$\begin{array}{c}\Delta r_{gr}=[\sqrt{1+{\gamma }_{pr}^{gr}B\tau }-1]r,\\ \\ {\gamma }_{pr}^{gr}=min(1,max(0,1-F_C/F_{cr}^{gr})),F_{cr}^{gr}=2f_S[r/r_0]{\sigma }_{cr},\end{array}$$

(13)

where $F_C$ is the sum of the projections of all forces affecting paired cells on the vectors, starting at the points of application of these forces and ending at the center of mass of the cell pair. Therefore, external pressure slows down tumor cell growth and a sufficiently large pressure can stop it. Growing cells need to consume nutrients at a rate proportional to their growth rate; therefore, for growing cells ${\gamma }_{pr}={\gamma }_{pr}^{gr}$.

If $n\le n_d$, cells die, which is described by the simple option, previously implemented in agent-based models (e.g., [32]), i.e., cells shrink, with the center of mass position remaining fixed until they disappear. If a cell radius $r>r_{min}$, it changes by the following negative amount:

$$\begin{array}{c}\Delta r_{sh}=[\sqrt{1-M\tau }-1]r,\end{array}$$

(14)

and when its radius reaches $r_{min}$, it disappears, meaning that it and its bonds and potential bonds are not accounted for from this moment on.

If the local nutrient level lies within the range $(n_d,n_q]$, tumor cells are considered to be quiescent and they do not shrink, grow, diverge, or die.

3.1.3. New Parameters

The additional parameters used in the simulations of encapsulated tumor growth are listed in

Table 2. Additional model parameters that were used in the simulations of tumor growth in tissue.

Parameter	Description	Model Value(s)
$k_N$	coefficient of stress function for normal cells	80/200/600
$k_C$	coefficient of stress function for bonds between cells and capsules	80
$K_C^T$	coefficient of tumor capsule speed of growth	5
$K_C^N$	coefficient of normal capsule speed of growth	0.5
$r_N^0$	initial radius of normal capsule	55/104/158/205
$s_0^N$	normalized stress needed to stop normal capsule shrinkage	0.1
B	tumor cells’ maximum rate of growth	0.03
M	tumor cells’ rate of shrinkage	0.01
$r_{min}$	radius of a tumor cell at which it dies	2.5
P	coefficient of nutrient inflow	1/4/16
$D_n$	nutrient diffusion coefficient	3000
$Q_N$	rate of nutrient consumption by a normal cell and a quiescent tumor cell	0.5
$Q_T$	maximum rate of nutrient consumption by a tumor cell due to proliferation	25
$n^{*}$	Michaelis–Menten parameter for nutrient consumption	0.005
$n_q$	nutrient level below which tumor cells become quiescent	0.1
$n_d$	nutrient level below which tumor cells shrink and die	0.01

. Note that the value of the tumor cells’ maximum rate of growth B was mentioned above when defining the tumor cell maximum divergence speed v.

The coefficient of stress function for the bonds between normal cells $k_N$ varied analogically to that for tumor cells. The corresponding value for the bonds between cells and capsules $k_C$ was chosen to be low in order to increase the free space for tumor cell proliferation. The parameters related to capsules were adjusted empirically in order to suppress excessive tumor outgrowth and tissue rupture but to inhibit tumor growth only rather slightly under a low number of normal cells. The tumor capsule is considered to be more compliant, as in reality, the growth of a tumor is associated with the remodeling of connective tissue due to the secretion of specific enzymes, which should make it more susceptible to deformation (see [47]). Experiments show that tumor growth can be significantly inhibited by the neutralization of these enzymes’ secretions [48]. The initial size of normal tissue was varied to study its effect on tumor growth. Under designated values of $r_N^0$, it contained 90, 329, 812, and 1380 normal cells. Tumor cells’ rate of shrinkage in relation to their death rate was roughly assessed based on the fact that for a lot of tumor cell lines, more than half of tumor cells can survive after two days of starvation [49].

The coefficients related to the nutrient were estimated based on the experimental data on glucose, which plays a crucial energetic role in tumor cell metabolism and is indispensable for cell biosynthesis [50]. These estimations were described in our previous works (see, e.g., [51,52]). The nutrient diffusion coefficient was, however, reduced by about 2.5 orders of magnitude compared to that of glucose to generate sufficient nutrient gradients within regions containing only hundreds and thousands of cells. The coefficient of nutrient inflow P has the physical meaning of microvasculature permeability surface area product. Both the permeability and surface area of microvasculature increase during tumor angiogenesis [53]; hence, the effects of angiogenesis and antiangiogenic therapy on tumor growth can be evaluated by varying this parameter.

3.2. Results

3.2.1. No Accounting for Nutrients

In this section, the simulations of tumor growth in normal tissue are discussed, without nutrient deficiency as an inhibitory factor. Tumor growth is affected by two types of stress—residual stress exerted by tumor cells on each other and external stress due to the displacement of normal cells, mediated by the tumor capsule.

Figure 6 shows the initial stages of encapsulated tumor growth within a small region of normal tissue, the initial radius of which is $r_N^0$ = 55$\mathsf{\mu }$m, for two values of tissue conductivity K. The upper row corresponds to a small K, resulting in a slow cell displacement in response to the applied forces. As a result, the approximate sixfold rotational symmetry of the normal tissue is maintained for a relatively long time during tumor growth. The reorganization of normal tissue begins with the displacement of six lines of normal cells, which are aligned perpendicular to the tumor capsule and originate from the six cells initially linked to it. These six cells, which form the first layer of normal cells around the tumor, diverge from each other in response to tumor growth and gradually integrate within the next layer of normal cells. All other cells within the six designated lines behave similarly. Their displacement is at first counteracted by the stretched bonds of the layers to which they initially belong, as the leftmost figure shows. Furthermore, their displacement is conversely facilitated by the newly formed stretched bonds with the cells of the next layers. This remodeling of a normal cell network results in a six-pointed structure, which begins to notably lose its symmetry only in the third week of tumor growth.

The lower row in Figure 6 corresponds to a high value of $K=10$, resulting in the cells displacing notably faster and local disturbances in the normal cell network being quickly redirected to remote areas. This results in much faster symmetry breaking due to unavoidable instabilities and leads to the formation of a more chaotic normal cell network, which transiently contains a notable number of void areas, as shown in the figure corresponding to 100 h of tumor growth. As shown in the figure on the right, during encapsulated tumor growth at $K=10$, a sufficient amount of cells in the outer layers divided at notably less than the maximum rate, contrary to the case of free tumor growth (see Figure 5), which suggests the slowing down of tumor growth by external stress.

Figure 7 illustrates the effect of the external stress on tumor growth by highlighting the times when the number of tumor cell pairs reached 3500 depending on the initial radius of the normal capsule $r_N^0$ for different values of the tissue conductivity K and the coefficient of stress function for normal cells $k_N$. The coefficient of stress function for normal cells $k_T$ is equal to 80 in the corresponding simulations so the interactions among tumor cells remained weak. The value of $r_N^0=0$ corresponds to the case of tumor growth in the presence of a tumor capsule but without normal cells so these points coincide for different values of $k_N$. Gray crosses on the ordinate axes denote the analogous moments for free tumor growth (corresponding to the results discussed in Section 2.2 minus $\approx 23$ h due to the alteration of the initial conditions).

At the highest value of K, which corresponds to the graph on the left, the presence of a tumor capsule in the absence of normal cells increased the time to reach 3500 tumor cell pairs by about $20\%$. At $k_N=80$ and 200, the presence of normal cells increased this parameter almost 2.5 times, and notably, the number of normal cells barely affected it within the considered range of $r_N^0$. This may seem counterintuitive at first since larger normal tissue should be more difficult to displace by a growing tumor. However, an increase in the number of cells also led to a more active formation of intercellular bonds, which emerged in the stretched state and tended to pull the normal cells into the less deep cell layers, as discussed above. This process conversely stimulated tumor growth. At the highest value of $k_N=600$, tumor growth was, however, halted at the very beginning under $r_N^0\ge 104$ since the stretched normal tissue exhibited a high enough pressure to stop tumor growth before tissue remodeling outbreaks.

At $K=0.32$, a more intuitive pattern of the decrease in the tumor growth rate with the increase in the normal tissue size manifested itself and became even more pronounced at $K=0.01$. This was due to the fact that the formation of new bonds between normal cells, which eventually stimulated normal tissue remodeling tumor growth according to the above-described mechanism, played a smaller role in the decrease in K (e.g., during the simulations depicted in Figure 6, ≈209 thousand different intercellular bonds were counted until there were 3500 tumor cell pairs at $K=10$ and only $\approx 92$ thousand at $K=0.01$). However, at $K=0.32$, tumor growth was halted from the very beginning with a significant size of normal tissue. This occurred as the cells and capsules became stuck in a stable position after the first several divisions of tumor cells. This effect can be attributed to the peculiarity of the model, which could be overcome by further alterations of the initial conditions; however, the corresponding modifications are not considered here since the results for extreme values of K are sufficient for obtaining the desired qualitative insights.

At $K=0.01$, the tumor grew up to 3500 cell pairs under all the considered values of $k_N$ and $r_N^0$. Higher values of $k_N$ resulted in faster tumor growth due to more active normal tissue remodeling. Contrary to the case of $K=10$, tumor growth did not stop at $K=0.01$ when normal cell interaction was the strongest, even under the greatest considered size of normal tissue. This was due to the much slower propagation of normal cell displacement along the tissue in this case, which allowed the tumor to grow, whereas the distribution of stress within the normal tissue had not yet stabilized. Conversely, at $K=10$, it stabilized very fast and deformations of remote stretched areas of normal tissue notably influenced the pressure exerted on the tumor capsule by the closest normal cells.

The results for higher values of $k_T$ are not shown, however, suffice it to say that the increase in this parameter significantly inhibited tumor growth and at $k_T=200$, tumor growth was halted at $r_N^0=104$ under any values of the other varied parameters. This was mainly due to the fact that under higher $k_T$, tumor cells have less available free space for division since the critical pressure is achieved at smaller distances between tumor cells. Therefore, overall, the model simulations suggest that the weakening of the interactions between tumor cells is a significant condition for tumor growth.

3.2.2. Accounting for Nutrients

In this section, nutrients are accounted for and nutrient deficiency influences tumor growth, along with solid stress. Figure 8 shows the dynamics of the normal capsule radius $r_N$ for the corresponding type of growth, with its initial value being $r_N^0=158$, for variations in the normal cell stress function coefficient $k_N$, tissue conductivity K, and coefficient of nutrient inflow P. The tumor cell stress function coefficient was $k_T=80$, as in the simulations described above.

At high tissue conductivity, $K=10$, tumor growth practically stopped during the first half a year under a moderate and a low nutrient inflow, $P=1$ and 4, and notably decelerated under $P=16$, whereas the variation in the normal cell stress function coefficient $k_N$ played only a small role in the outcome. The tumor, which was stable in size, is shown in the bottom-left picture, along with the normal tissue that was torn in two from the pressure exerted by the growing tumor. As nutrients are supplied from normal tissue, the configuration of the normal tissue also affects the tumor. Therefore, the necrotic tumor core, devoid of cells, has an elongated shape and is outlined by dying cells, which are depicted in lighter shades. As these cells shrink and are eliminated, they pull other tumor cells into the tumor core. The ongoing proliferation of tumor cells in the outer layers, where the nutrient concentration is sufficient, compensates for the death of those cells. Within such a dynamic balance, the number of tumor cell pairs fluctuates at around 800 under the corresponding values of the model parameters.

At a low tissue conductivity $K=0.01$, tumor growth was rapidly halted at $P=1$. However, it continued for a long time, at least 5 years, under sufficiently high nutrient inflow and high $k_N$. The latter, as discussed above, yielded rapid remodeling of normal tissue in response to its perturbations. An example of the corresponding tumor is shown in the bottom-right picture. At the depicted moment, the normal tissue had already stretched enough to be only one cell in width, and the tumor core already occupied the bigger part of the tumor volume. However, the nutrient inflow was sufficient to keep ≈8000 tumor cell pairs alive, their number slowly growing as the tumor radius increased. With that, the stress generated by dying tumor cells was insufficient to notably slow down the tumor growth at such a low tissue conductivity.

4. Conclusions

This work introduced an off-lattice agent-based model of tumor growth, where the tumor is represented by an interconnected network of proliferating cells. The rate of proliferation of each cell depends on the solid stress generated due to its interaction with other cells. Sufficiently closely located cells tend to repel from each other, whereas with the increase in the intercellular distance, this interaction at first weakens and then turns to adhesion, which strengthens up to a certain intercellular distance and then also weakens. An original feature of the model is the consideration of tumor cells as pairs of overlapped diverging circles, with a specific method used to account for intercellular bonds, which ensures smooth dynamics of the cell network and maintains relative numerical efficiency.

Two different types of tumor growths were considered during the simulations. The first was free tumor growth, restricted only by residual stress generated within the tumor. The second was tumor growth within a round capsule surrounded by normal tissue, represented by interconnected non-proliferating normal cells, surrounded by another capsule. The capsules were introduced to simulate the elastic response of a tissue, which, in reality, is provided by extracellular matrix elements. The tumor growth in tissue was simulated both with and without accounting for nutrients supplied from the normal tissue that influence the behavior of tumor cells.

The simulations of these growth regimes led to the following conclusions:

• For free tumor growth, it was demonstrated that tissue hydraulic conductivity, which determines the speed of cell displacement in response to stress, is the major factor that defines the rate of tumor growth. Tumor growth can be close to exponential only at sufficiently large values, whereas under low values, the cells within the tumor volume experience significant stress that stops their proliferation. The form of the freely growing tumor crucially depends on the strength of cell–cell interaction, as strongly interacting cells tend to cluster, yielding an infiltrative tumor shape.
• When encapsulated tumor growth in normal tissue was considered, the strong interaction of tumor cells turned out to be a major factor that significantly limited tumor growth even within a thin normal tissue. This is in agreement with the experimental data that showed that the process of carcinogenesis leads to the weakening of the strength of intercellular contacts [54]. The future development of the model will include the consideration of heterogeneous tumor cell adhesive properties and their evolution due to the random alterations upon cell division, which can be quite naturally reproduced through an agent-based approach.
• At high tissue conductivity, another major factor that led to the rapid halting of the growth of an encapsulated tumor in normal tissue was the strength of the interaction between normal cells. High values led to the inability of the growing tumor to initiate normal tissue remodeling, which is necessary to provide the space for tumor growth. However, the simulations at low tissue conductivity showed that once remodeling of normal tissue was initiated, the strong interaction between normal cells nevertheless resulted in the acceleration of this process, stimulating fast tumor growth.
• An increase in the initial size of normal tissue at low tissue conductivity led to a decrease in the tumor growth rate as expected. However, at high tissue conductivity, the increase in the number of normal cells barely influenced the tumor growth rate due to more active remodeling of remote areas of normal tissue. This observation was made for small sizes of normal tissue of no more than $0.4$ mm in radius, whereas the verification of this phenomenon for larger tissue sizes requires significant computing power. To avoid high computational costs, a continuous model could be created whose dynamics correspond to those of the presented agent-based model and could be used to qualitatively reproduce its results. This lies within the scope of future work.
• It is worth noting that the results obtained when considering tumor growth in normal tissue when accounting for nutrients are in good qualitative agreement with the previous results obtained using a simple continuous model in which tumor and normal tissue behaved as an elastic fluid-like substance [52], where sufficiently low tissue conductivity and sufficiently high nutrient supply in both models led to continuous tumor growth for at least several years, with the formation of large necrotic cores. Notably, the formation of giant long-growing benign tumors is indeed a clinically observed phenomenon in cases of tumors arising from connective tissue with inherent low tissue conductivity [55,56].