The tumor as an organ: comprehensive spatial and temporal modeling of the tumor and its microenvironment

Background Research related to cancer is vast, and continues in earnest in many directions. Due to the complexity of cancer, a better understanding of tumor growth dynamics can be gleaned from a dynamic computational model. We present a comprehensive, fully executable, spatial and temporal 3D computational model of the development of a cancerous tumor together with its environment. Results The model was created using Statecharts, which were then connected to an interactive animation front-end that we developed especially for this work, making it possible to visualize on the fly the on-going events of the system’s execution, as well as the effect of various input parameters. We were thus able to gain a better understanding of, e.g., how different amounts or thresholds of oxygen and VEGF (vascular endothelial growth factor) affect the progression of the tumor. We found that the tumor has a critical turning point, where it either dies or recovers. If minimum conditions are met at that time, it eventually develops into a full, active, growing tumor, regardless of the actual amount; otherwise it dies. Conclusions This brings us to the conclusion that the tumor is in fact a very robust system: changing initial values of VEGF and oxygen can increase the time it takes to become fully developed, but will not necessarily completely eliminate it. Electronic supplementary material The online version of this article (doi:10.1186/s12859-016-1168-5) contains supplementary material, which is available to authorized users.

The Area object is the controller of the system. It describes the overall "world" of the system and handles the three dimensional lattice of the objects' positions within the world, the molecules' positions, and the time steps of the system, as well as taking care of outputting the data to files. It contains many functions that are used throughout the simulation.
Each of the objects is also linked to a Position entity, which holds their specific and current x, y, z coordinates.
When the model is executed, the dynamics of the system emerge as a product of the behavior of all the entities and the interactions between them.

Cell size
In order to create a realistic tumor, there was a need to properly define the position occupied by a cell. If it were to take up a single position the tumor would take on a symmetric form. To overcome this problem, we defined the resolution of the cell to be a cube of any edge length in terms of the number of positions, where if it is above 1, the tumor will take on a non-symmetric form. This is due to the fact that a neighboring cell does not necessarily have to be an exact position above or below it, but can also be half a cell-size above (in the case of cell resolution of 2), etc.

Tumor cell survival & proliferation
The tumor cells can only proliferate under the condition that there is free space, the size of a cell, around them. The position of the new cell is chosen randomly from any adjacent place around the cell (including the diagonals). A tumor cell constantly consumes oxygen and will proliferate only if it has a sufficient amount of oxygen.
The cell consumes the oxygen it needs randomly from the cell and its immediate surroundings of one position in any direction around it. If it does not have a sufficient amount, it will become hypoxic (deprivation of adequate oxygen supply) and stop proliferating. In parallel it will start sending angiogenic factors, VEGF. If it continues to lack sufficient amounts of oxygen it will become anoxic (extreme decrease in the level of oxygen) and eventually die (become necrotic). An active hypoxic cell that returns to consume enough oxygen will stop secreting VEGF.

VEGF diffusion
VEGF diffusion is a continuous process, but we chose to represent it in its most basic formas individual molecules that move in a random walk fashion. Each molecule that is secreted takes on the initial location of the tumor cell from which it came, and takes up exactly one position. It can move into any adjacent position (including all diagonals), even if that position is taken up by another molecule or cell. Each VEGF object can also represent more than a single molecule. While implementing this, there was a need to check if the collection of the individual molecules moving randomly indeed creates the overall diffusion phenomenon. This involved testing the amount of molecules that each cell secretes, the rate of secretion, and the rate of the molecules' random walk, and also deciding on what happens to molecules that reach the edge; should they disappear, come back, or stay at the edge? These tests resulted in a successful diffusion progression.
Every VEGF molecule in the model represents ~10 real VEGF molecules: In [1] VEGF values in tumor were measured. From table 6 in this paper, the total tumor VEGF that was measured was 13 pg/mg tumor. This value was converted to 0.3 pmol/cm 3 tissue to be used in their model. Since the world of our model is of size 100*100*100 positions, where each position represents 5 µm, it is 500 3 µm 3 = 0.05 3 cm 3 = 0.000125 cm 3 . Therefore there will be 0.3*0.000125 pmol in our model's world = 0.0000375pmol = 3.75 *10 -5 *10 -12 mol = 3.75*10 -17 *6 *10 23 mol = 22.5 million molecules. In our model, at the stage when the tumor is developed and vascularized the average number of VEGF molecules is 2-3 million, meaning that each molecule in the model represents ~10 real VEGF molecules. This is a reasonable number, allowing, on the one hand to represent only a fraction of the real number, but only to one order of magnitude, meaning that it is sensitive enough to observe changes that occur to the molecules. Every oxygen molecule in the model represents ~7 million real oxygen molecules: In [2], 10mmHg was found to be the median partial pressure of O2 in tumors. By using the ideal gas equation of PV=nRt, where P=pressure in pascals = 10mmHg*133 = 1333pascals, V=volume in cubic meteres, in our model, the world is 500 3 µm 3 = 0.0005 3 m 3 = 0.000000000125 m 3 = 1.25*10 -10 m 3 , R (constant) =8.314, T=temperature in kelvin = 37 degrees Celsius in the body = 310 Kelvin. Therefore n (number of mols) = 1333*1.25*10 -10 / (8.314*310) = 6.46*10 -11 mols = 6.46*10 -11 *6*10 23 =~39*10 12 oxygen molecules in the system. This corresponds to ~5-6 million molecules on average in the model, meaning that each oxygen molecule in the model represents ~7 million molecules. This too is a reasonable number, as there are many more oxygen molecules, as they are much smaller and are largely consumed by all cells everywhere all the time for survival. It is obviously impossible to represent them individually and their gradient is less important than with the VEGF molecule, but more importantly is to be able to observe where enough oxygen is present.

Molecule elimination
Both VEGF and oxygen that pass the edge of the world or are consumed by the appropriate cells, are eliminated from the system.

Angiogenesis
The angiogenesis process in the model occurs along the following steps:  Endothelial cell checks for VEGF around it (environment sensing).
 If VEGF is found it binds to the cell according to a binding probability, and the molecule is eliminated.
 If the levels of VEGF molecules that bind to the cell pass a certain threshold (the angiogenic switch), the cell becomes activated and begins the process of angiogenesis.
 This process is guided by the VEGF gradient that the cell senses at each point.

Angiogenic switch
In order for the endothelial cell to be activated for angiogenesis it needs to pass the threshold of meeting a certain amount of VEGF's in a certain amount of time [3].

Endothelial direction
The direction taken by the endothelial cells during angiogenesis has a great effect on the final spatial organization of the vessels. In our model, the cell scans a vector radius around it and calculates the geometric averages of the amount of VEGFs and moves in the direction of the highest average. The radius scanned can be increased or decreased to change the manner of the vessel's elongation; the larger the radius being scanned the more directly the vessel will progress towards the tumor, and the smaller the radius the more twists and turns it will have while progressing towards the tumor in an indirect manner.

Endothelial elongation
If the endothelial cell is directed by the VEGF gradient to elongate into a position occupied by another endothelial cell, it will join this vessel and stop elongating. If, however, the position is occupied by a tumor cell, the endothelial cell will try to find an empty space in the closest possible direction, so as to continue its elongation. In addition any endothelial cell that is a part of an angiogenic vessel can also split and branch out of its main vessel if it binds to a high amount of VEGF in a specified time.
The continuation of the endothelial cell progression depends on the presence of VEGF molecules. Without a sufficient amount thereof it will not continue to proliferate and the vessel will halt the elongation.

Delta-Notch inhibition
Delta-notch lateral inhibition signaling is involved in cell-to-cell communication and regulates the determination of various cell fates. Among cells that have the potential to adopt the same fate, lateral inhibition specifies some cells for a primary or preferred fate and others for a secondary or alternative fate [4,5]. In the case of endothelial cells, the activated cell inhibits its adjacent cells, thus preventing them from becoming active. A list of the parameters used in the cancer model together with the description of each parameter and its default value. The default values are used to achieve a standard run and can change according to changes in other parameters, as well as due to further development of the model. Parameters are unitless but are quantified relative to each other. Each time step (ts) in the model is approximately equal to 1-2 hours, and each 10 3 cells in the model are approximately 100 3 actual cells.