Mathematical simulation and prediction of tumor volume using RBF artificial neural network at different circumstances in the tumor microenvironment

Uncontrolled proliferation of cells in a tissue caused by genetic mutations inside a cell is referred to as a tumor. A tumor which grows rapidly encounters a barrier when it grows to a certain size in presence of preexisting vasculature. This is the time when it has to find a way to go on the growth. The tumor starts to secrete tumor angiogenic factors (TAFs) and stimulate preexisting vessels to grow new sprouts. These new sprouts will find their way to the tumor in the extracellular matrix (ECM) by the gradient of TAF. As these new capillaries anastomose and reach tumor, fresh oxygen is available for the tumor and it will reinitiate the growth. Number of initial sprouts, distance of initial tumor cells from the vessel(s) and initial density of the tumor at the time of sprout formation are questions which are to be investigated. In the present study, the aim is to find the response of tumor cells and vessels to the reciprocal effects of each other in different circumstances in the tissue. Together with a mathematical formulation, a radial basis function (RBF) neural network is established to predict the number of tumor cells at different circumstances including size and distance of initial tumors from the parent vessel. A final formulation is given for the final number of tumor cells as a function of initial tumor size and distance between a parent vessel and a tumor. Results of this simulation demonstrate that, increasing the distance between a tumor and a parent vessel decreases the number of final tumor cells. Specially, this decrement becomes faster beyond a certain distance. Moreover, initial tumors in bigger domains must become much bigger before inducing angiogenesis which makes it harder for them to survive.


Introduction
Cancer is known as the second leading cause of death worldwide. 1 Therefore, understanding the mechanisms of this uncontrollable disease is of vital importance. Rapid proliferation of these abnormal cells depends on different spatial and temporal scales 2 that is affected by the tumor microenvironment and the structure of the surrounding tissue. 3 As known, tumor growth includes three distinct levels. 2 At first, a tumor grows while it is fed by existing vasculature in the host tissue. Since tumor cells proliferate rapidly and consume more oxygen, the surrounding vasculature is not able to meet nutritious requirements of the starving tumor cells. In the second stage, the tumor seeks a way to overcome oxygen deficiency by inducing vessel growth toward the tumor. The process of stimulating and growing new capillaries from pre-existing vessels is known as angiogenesis. The third and most fatal part of the growth is when the tumor has acquired its own vascular network. Unstable and leaky capillaries are not only a source of oxygen, but also a shortcut for tumor cells to breach into and be transported to the other parts of the body which is called metastasis. Understanding these distinct levels and dependence of tumor on surrounding tissue is crucial in defining the outcomes of a chemotherapeutic treatment.
For many years, scientists were struggling to find the underlying mechanisms of this devastating disease. 4,5 Folkman was the first scientist who distinguished the different steps of tumor growth. 6 He then proposed the hypothesis that tumors are not able to grow beyond a certain size in the absence of angiogenesis. 7 Since then, many scientists and clinicians have been working on angiogenesis and vascularized tumors to fully capture the underlying phenomena and hinder tumor growth via obstructing the ways for tumor to provoke angiogenesis. A momentous part which is still obscure for scientists is the effects of surrounding tissue and existing vasculature. In vivo models are difficult and costly to implement. Moreover, all physical phenomena happening in the tumor microenvironment cannot be investigated by in vivo models. On the other hand, in vitro models make clinicians able maneuver on some precise and controllable parameters. For instance, it is possible to culture endothelial cells and to see if they respond to some stimulating factors or not. However, these models have drawbacks which make them limited in a variety of ways. Restrictions related to temporal and spatial scales inhibit utilization of experimental methods for tumor studies. 8 The only choice remaining is using in silico models which make researchers able to test their hypotheses without considerable effort in comparison to experimental methods. In addition, when new findings are acquired in experiments, scientists refer to in silico models to clarify and illuminate the underlying stimulators for the phenomena. Therefore, scientists and researchers have recently been attracted toward mathematical models and it is believed that integrated methods are the most reliable and efficient ways to study tumor growth. 9 Wide variety of methods, parameters, physical situations, and time scales could be modeled through in silico models.
According to Araujo and McElwain, 9 Hill 10 was pioneer in studying the tumor microenvironment mathematically. Following him, many scientists have studied avascular and vascular tumor growth using different approaches leading to important outcomes which have assisted clinicians in the battle against their nemesis. All in all, mathematical models are categorized into three main groups of continuum, discrete, and hybrid. Each method is accompanied by its pros and cons. Continuum methods are governed by Partial Differential Equations (PDEs) and treat tumor as ensembles of cells acting and moving together. Hence, this modality is able to capture macroscale interactions between tumor cells. On the other side of mathematical models lies discrete ones. These methods as known by their names, handle each single cell inside the tumor microenvironment. As a result, it is feasible to track each cell and decide about its destiny in the model. Nonetheless, in this case, some macroscale phenomena cannot be simulated. Moreover, as the number of cells increases, time, and cost of simulation increases concurrently. To overcome these difficulties, researchers have combined two methods and hybrid formulation is achieved. In this formulation, some part of the simulation is governed by PDEs and the other by discrete methods. In some cases, it is suitable to use continuum methods to capture the phenomena in macroscale which is costly and difficult by discrete models.
Different macroscale aspects of tumor growth have been investigated by continuum methods. 1,[11][12][13][14][15] Mechanical interactions, cell-cell adhesion, cell-matrix interactions, and ECM stiffness are studied by the so called method which occur in macroscales. On the other hand, most models use discrete framework to treat tumor growth numerically. 2,[16][17][18][19] These methods have reported intracellular and intercellular phenomena such as cell birth and death, cell deformation and cell-cell signaling. And also hybrid models which are a combination of both models is being used to study tumor. Scientific literature is rich of papers about different mechanical and chemical interactions among cells and the tissue. Dyson et al. 20 studied the mechanical properties of ECM and cell-ECM interactions and fiber anisotropy on the morphology of tumor cells. Authors stated that different patterns for cell aggregates and collagen and medium were observed under different mechanical conditions. The results showed that fibers realign by cell migration. In addition, circular, elliptic, and stripe configurations for tumor cell aggregates were observed. Cho and Levy 11 evaluated the treatment efficacy of cytotoxic and cytostatic drugs. They concluded that tumor cells were totally eliminated when both drugs were infused in the tissue simultaneously and relapse is prevented. Salavati et al. 21 worked on a new formulation of vascular tumors and validated the results by experiments on mice. The model confirmed the necessity of angiogenesis in solid tumors and the model can be used as a basis formulation to study therapeutics in solid tumors. Santagiuliana et al. 22 assessed the biomechanical effects of the tumor microenvironment on tumor growth over time. According to the authors, decreasing ECM stiffness and cell-ECM adhesion and increasing porosity of ECM increases tumor mass. Sefidgar et al. 23 studied effects of ECM stiffness on Endothelial Cell (EC) migration speed and reported that the ECs migrate faster in intermediate ECM stiffness. Salavati and Soltani 24 investigated impacts of proliferation rate on the overall rate of ECs movement using a Cellular Potts Model (CPM). They concluded that using variable proliferation rate for ECs, more realistic results could be achieved which matches well with experimental observations. Effects of tumor shape and size on drug delivery is studied by Soltani and Chen. 25 The results where indicative of importance of tumor size in filtration flux (flow rate out of the vasculature per unit volume) and pressure inside tumor. Bazmara et al. 26 studied EC proliferation and migration focusing on the effects of blood flow on EC phenotype alteration. The results showed that in the absence of blood flow and shear stress, a loop cannot maintain its stability and collapses at the end. Following this research, Bazmara et al. 27 published a paper studying blocking of intracellular signaling pathways to inhibit lumen formation and EC proliferation and hence to prevent angiogenesis during tumor growth.
Smeared finite element method which is introduced by Kojic et al. 28,29 and further studied by Milosevic et al. 30 can be used to model mass transport between capillaries and tissue in the tumor microenvironment. Kremheller et al. 31 used smeared representation of neovasculature to capture species transport in tumorinduced vasculature which is highly tortuous and contains tiny vessels. They compared the results obtained by their model to preexisting models to validate their results.
Computational studies have been beneficial to assess the efficiency of chemotherapy for cancer treatment. Ribba et al. 32 used an age-structured cell cycle model to examine efficacy of anti-metastatic agents, called inhibitors of matrix metalloproteinases (MMPi). It was revealed that MMPi is not an effective treatment in advanced cancers. Billy et al. 33 coupled a continuum model of angiogenesis with a continuum model of tumor growth to develop an age-structured cell cycle model and examine influence of endostatin overproduction on hindering angiogenesis. It was concluded that endostatin overproduction may suppress angiogenesis and tumor growth. Ribba et al. 34 investigated the effect of chemotherapy on tumor growth using an age-structured model and demonstrated that physiopathological parameters have higher impact on treatment efficacy than drug-related parameters. Lignet et al. 35 studied the interaction of antiangiogenic drugs and chemotherapeutic agents and presented an optimal therapeutic strategy for combined usage of these two treatments.
Pamuk et al. 36 proposed a mathematical model to study tumor growth. They studied onset of vascularization of tumors inside the tissue. Their results were in good agreement with experiments. Perfahl et al. 37 studied effects of domain size and boundary conditions on tumor response. Their results proved that distance of tumor and parent vessel, initial structure of the vascular network are important in determining the final architecture of the vascular tumor. Grogan et al. 38 performed a mathematical model of mouse cornea to study effects of cornea geometry on vascular network formation. They studied effects of geometry on VEGF concentration and vascular density in the domain. The results proved that selection of pellet location and geometry of the domain greatly influences angiogenesis. Modeling of avascular tumor growth, blood flow and vessel regression are absent in their modeling. Recently, Nikmaneshi et al. 39 proposed an agent-based mathematical model to study avascular and vascular tumors. The authors surveyed abnormal blood conditions encountered in cancer patients. It is concluded that hyperglycemia, hyperoxemia, and hypercarbia can aid in tumor growth. Moreover, hypertension can also lead to malignant tumors.
One important feature of tumor growth which is a debate among scientists in spite of great advancements in in silico models is the behavior of tumor cells in different distances of parent vessels. Which environmental conditions can make cells more devastating and when clinicians can better suppress tumor growth during the growth? Gimbrone et al. 40 experimentally studied tumor growth in different conditions of parent vessel. Their results were helpful in understanding the states and reaction of tumor cells in the tissue. Recently, Ghazani et al. 41 studied growth of a solid tumor at different distances of a parent vessel using a discrete tumor model where effects of tumor interstitial pressure in tumor cell cycle are also included.
In this paper, the aim is to make a useful contribution in understanding the fatality of a tumor in a variety of situations using a continuum tumor model and artificial neural network.
Effects of tumor location in the domain, domain size, and number of initial sprouts on the final architecture of the tumor and capillary network using a continuum tumor model and an agent based probabilistic angiogenesis method are studied.

Mathematical model
In order to model growth of a tumor in a tissue, a hybrid continuum-discrete model of tumor growth is prepared. The mathematical model for simulation of tumor cells is a continuum age structured model which makes use of PDEs. A discrete agent-based method which comprises probabilistic phenomena is used to make us able to track each single EC in the domain. Therefore, branching, anastomosis and blood flow are considered in the model. Based on blood flow, oxygen diffuses in the domain and reaches nearby tumor cells. On the other hand, VEGF as dominant TAF is secreted by tumor cells and consumed and used by ECs to route their way.

Tumor growth
Tumor growth is modeled by a continuum agestructured cell cycle model embedded in a macroscopic model of tumor dynamics. 33 Four different cell types are in tumor microenvironment: Proliferative, quiescent, necrotic, and healthy cells. To model mitosis, proliferative cells are considered to be in two different phases named proliferative phase 1 and proliferative phase 2. In ''Restriction Point,'' oxygen concentration and available spaces are evaluated to define where tumor cells are placed in their cell cycle. Cells enter proliferative phase 2 when oxygen concentration is sufficient for proliferation and space is available for cells to duplicate. As soon as cells pass the restriction point and enter proliferative phase 2, they continue their cell cycle disregarding oxygen availability. If defined conditions are not satisfied, cells are transferred to quiescent state. Cells can survive unlimited in the quiescent state until the conditions are met. If oxygen concentration falls below severe hypoxia threshold, the cells cannot endure severe oxygen deficiency and become necrotic (the condition is modeled as function h). Therefore, proliferative cells develop in both age and time in their cell cycle while quiescent and necrotic cells develop only in time. Figure 1 shows schematic of cell cycle.
To model transition between different phases in tumor microenvironment, functions g and h are introduced as follows in point (x, y) and time t in equations (1) and (2).
where P 1 , P 2 , Q are proliferative phase 1, proliferative phase 2 and quiescent cell's density, O is oxygen concentration, t O , t 1,h, and t 2,h are overpopulation threshold, hypoxic threshold and severe hypoxic threshold, a is the cell's age in time, and a max , P 1 , a max , P 2 are maximum duration of proliferative phases P 1 and P 2 , respectively. Since there are two proliferative phases in the model, total density of proliferative phase is calculated in point (x, y) and time t in equation (3).
where P is the total density of cells in proliferative phase. Since the total number of cells per unit volume is constant, densities are non-dimensionalized by constant total number of cells, N 0 , as shown in equation (4).
where A and M are densities of necrotic and healthy phases, respectively. Tilde, (;), is omitted in equations (5) and (6) for the sake of brevity. Tumor cells develop in both age and time in proliferative phases, therefore the evolution equations are presented in equations (5) and (6).
Quiescent and necrotic cells develop in time only in equations (7) and (8).
where g 0 t ð Þ + , h 0 t ð Þ + and g 0 t ð Þ À , h 0 t ð Þ À are positive and negative parts of the derivatives, and v P 1 , v P 2 , v Q , v A are corresponding velocities of proliferative phase 1, proliferative phase 2, quiescent phase and necrotic phase, respectively. Q t À ð Þ implies the density of quiescent phase in previous time step. Boundary conditions for turning one phase to another are shown in equation (9).
The cell birth provides a passive movement and activates tumor cell motility. Given the fact that cells move due to passive movement created by the cell birth and are pushed out through the domain boundaries, equation (10) is dedicated to evolution of healthy cells.
The duplication of cells increases pressure locally and forces cells to move and relax their pressure. The Darcy law for porous media is applied to relate the pressure created by the cell birth and the velocity of tumor cells in tumor microenvironment (equation (11)).

v = À ru ð11Þ
in which v is velocity of tumor cells and u is a potential corresponding to the pressure in tumor microenvironment. It is assumed that all tumor cells move with the same velocity as v Combining equations (5)- (8) and (10) provides the relationship shown in equation (12), which is an indicator of tumor expansion.

Angiogenesis
Tumor-induced angiogenesis, growth of new capillaries from pre-existing vessels, was almost neglected until the time when Folkman 7 proposed the idea that angiogenesis is the most fatal part of a tumor's growth. Many scientists have tried to model angiogenesis based on mathematical models. The one which was proposed by Anderson and Chaplain 42 has attracted more attention and used by many researchers to simulate angiogenesis. The angiogenesis model in this paper also makes use of formulations represented by Anderson and Chaplain. 42 Based on Anderson and Chaplain, 42 there are three mechanisms governing ECs movement in the domain. Random movement which is a consequence of random walk of ECs which has led to a formulation similar to diffusion mechanism. The other mechanism affecting ECs in the extracellular matrix (ECM) which is proved by Keller and Segel 43 is chemotaxis, interpreted as direct movement of ECs as a result of soluble chemoattractants. The remaining mechanism in this model is haptotaxis, which is representative of transverse movements and is a result of adhesion force in the ECM. Combining these mechanisms and related formulations, total flux term for EC density in the domain is reached as equation (13): In which J n is total flux of ECs, J random , J chemo , and J hapto are diffusion, chemotaxis, and haptotaxis terms, respectively. Simplifying and writing the related parameters for each term results in equation (14) as below: where n is density of ECs, D n is diffusion coefficient of ECs, x(c) is chemotactic function, c is VEGF concentration, r 0 is positive coefficient for haptotaxis and f is fibronectin concentration. To nondimensionalize equation (14), the related parameters are divided by reference parameters introduced in Anderson and Chaplain 42 asc = c c 0 ,f = f f 0 ,ñ = n n 0 , andt = t t . Using these nondimensional parameters, equation (14) is nondimensionalized and a new dimensionless equation is reached: Nondimensional parameters in equation (15) are: Most of the parameters are adopted from Anderson and Chaplain 42 except the parameters which are related to domain size and should be updated to fit this model. t is dimensionless time parameter which is dependent upon the distance between the tumor and parent vessel and should be varied based on tumor position in the domain. The related dimensionless time parameters for different distances are stated in the Supplemental File. Discretizing equation (15) using Euler finite difference scheme, equation (17) is reached: An EC's movement is decided based on this equation in which R 0 is the probability to stay stationary and R 1 , R 2 , R 3 , R 4 are probabilities to move right, left, up and down, respectively. l, m are indicators of position of a point in the discretized computational grid in the direction of x and y, respectively. These probabilities are normalized by the sum of the probabilities. Then each movement probability falls in a range between 0 and 1. Therefore each limited range represents a movement to one of the four directions and also staying stationary. A random number is generated between 0 and 1 to decide about the direction of movement. Therefore, EC moves in the direction where the random number has fallen.
As stated above, fibronectin exists in the ECM. It adheres to the matrix and does not diffuse in the tissue. Therefore, its concentration is constant in the tissue initially. It is proved that fibronectin is secreted by ECs and also exists in the domain. On the other hand, fibronectin is consumed as ECs move in the domain toward the tumor. Summing all these factors, equation (18) is reached in which there are generation and consumption terms in the formulation: where v and m are positive constants. For ECs to move in the domain, a gradient of VEGF is required. The ECs sense the gradient and move toward the highest gradients. VEGF is consumed while ECs move in the domain which affects their subsequent move. Therefore, the equation of consumption of VEGF is considered in angiogenesis section as below: in which l is a positive constant. Secretion of VEGF by tumor cells is proposed in the next sections.
Equations (18) and (19) are nondimensionalized using the aforementioned parameters of reference EC density, VEGF, fibronectin, and time constants: where h, b, and g are positive non-dimensional coefficients defined as: The set of equations (20) is also discretized using Euler finite difference method.
Posterior to cell movement, ECs may encounter each other during movement and fuse together, which will lead to anastomosis. An EC can also come across a stalk cell and anastomosis will occur. On the other hand, ECs may branch in their way toward tumor. Many factors affect EC branching which can be introduced as: (1) ECs and stalk cells must reach a certain age and be mature enough prior to being able to branch, (2) Density of ECs in the point of branching is high enough, and (3) there is an available space nearby for the new EC to occupy. If these conditions are satisfied, one another parameter will be checked to let ECs branch, concentration of TAF must be above a threshold value to let ECs branch.
The probability condition is estimated in this paper to make the model able to mimic real physical situations in the tumor microenvironment. As stated above, ECs encounter each other and form closed loops. Circulating blood flow is the main aim of tumorinduced angiogenesis to deliver oxygen and nutrients to the tumor site. Therefore, as closed loops form, blood flows in the network and provides fresh oxygen for tumor cells. Dewhirst et al. 44,45 stated that there is no relationship between vessel diameter and blood flow inside tumor-induced neovessels. Nevertheless, many scientists 39,46-48 have assumed that since blood flow is laminar inside the neovessels, one can use Poiseuille's law from Mechanical Engineering formulations to obtain the flow in the vascular network. Therefore, in this research the same procedure is followed and blood is assumed to flow laminarly in the capillaries and obey Poiseuille's law: where DP B is pressure difference between two junctions, R is vessel radius, L c is distance between two junctions, and m B is blood viscosity. The blood viscosity is assumed constant and equal to 4 3 10 23 Pa.s. 49 Apart from blood flow rate, pressure difference is also unknown in this equation, hence another equation is required to determine the unknowns. Conservation of mass law is applied for each junction of the vascular network where three capillaries meet.
where N is the number of capillaries joining junction in the position shown by indices (l, m). Pressure difference across parent vessel is 8000 Pa (60 mmHg) and diameter of parent vessel is 14 mm and new capillaries are 8 mm. 49

Vascular growth
Two distinct parts of tumor growth are dependent upon each other via communicating through chemicals in the tumor microenvironment. Among many factors being transmitted via ECM between tumor cells and vessels are nutrients and TAFs. The most dominant TAF here is Vascular Endothelial Growth Factor (VEGF) which is secreted by hypoxic tumor cells and consumed by ECs. On the other hand, vessels are sources of oxygen which is considered to be dominant nutrient in the model. As new capillaries anastomose and blood circulates in the domain, the required oxygen for hypoxic tumor cells are provided and the cells will enter proliferative phase if the conditions are suitable. Flowchart of simulation of vascular tumor growth is demonstrated in Figure 2.
Dissemination of these chemicals are so fast in comparison to tumor growth and angiogenesis phenomena that the equations for concentration of them can be considered to be in steady state conditions. For distribution of VEGF in the domain, equation (25) is proposed as below: where D V is VEGF diffusion rate in the tissue, a V is production rate of VEGF by quiescent cells, and d V is degradation rate of VEGF in the tissue. There is no consumption term in this equation since it is considered in equation (19) where ECs consume VEGF when they move in the domain. The governing equation for calculating oxygen concentration is defined as equation (26).
,O are consumption rates of proliferative, quiescent, and healthy cells, respectively. d O is degradation rate of oxygen within the tissue domain. Required parameters are tabulated in Table 1.

Numerical method
To explore the effects of domain size (corresponding to tissue size), two different sizes of the domain are selected. The small domain is selected to be 4 mm 3 4 mm with 200 3 200 grids and the large domain is selected to be 6 mm 3 6 mm with 300 3 300 grids. For ECs with the diameter of 8-12 mm, the average diameter of 10 mm has been selected.
The discrete part of equations that belongs to angiogenesis are discretized using Euler finite difference method. The configuration of domains provides the mesh size of Dx = Dy = b = 20 mm. ECs are assumed to move one double space in real tissue during each time step of the simulation to enable us to use the same computational grid for the equations of tumor growth and angiogenesis. Time step is set equal to 1 h when we multiply time step to non-dimensional time parameter. A meshing network of 300 3 300 is used to preserve the size of meshes for the large domain. Therefore, every procedure explained for the small domain is also valid for the large domain.
No-flux boundary condition is imposed on ECs for angiogenesis equation (equation (27)).
where z is an outward unit normal vector. Finite volume method is used to solve equations of tumor growth. Cell-centered method is used for discretization of equations. The small and large domains consist of 200 3 200 and 300 3 300 control volumes, respectively. Control volumes are chosen such that center points of control volumes coincide with the points used to solve discrete equations of angiogenesis and eliminate the need to interpolation among the points. Since there is VEGF consumption rate (20) no specific restrains for selecting time step of the finite volume method, for the distances of 1, 1.5, and 2 mm time step is equal to 1 h and for the distances of 2.5 and 3 mm, time steps are 2 and 3 h, respectively, when we multiply the nondimensional time step to the nondimensional time parameter. Hence, in order to solve the coupled equations, for each step of tumor cell growth in cases of 2.5 and 3 mm distance, we go through 2 and 3 steps in angiogenesis, respectively. Boundary conditions for transition among phases are detailed in section 2.1. For the outer boundary, it is assumed that healthy cells leave the domain as tumor grows.
Oxygen and VEGF are used to couple two distinct parts of tumor growth. Steady state equations are used for oxygen and VEGF since their diffusion time scale is much shorter than the time scale of tumor growth and angiogenesis. Equations (25) and (26) are solved using Gauss-Seidel iteration method in the domain. We assume that oxygen and VEGF cannot cross the boundary, therefore Neumann boundary condition is imposed on the edges of the domain for oxygen and VEGF ( ∂O ∂z j ∂O = 0 and ∂c ∂z j ∂O = 0). Each configuration has simulated at least for five times and the mean values are reported in the results section.

Initial conditions
Tumors with a quiescent core and a proliferative rim are considered as initial tumors for simulations. In proliferating rim, outer strip is in proliferative phase 2 and inner strip is in proliferative phase 1. To realize how the size of initial tumor affects angiogenesis process and tumor growth is influenced by angiogenesis, tumor growth and angiogenesis are studied for three initial tumor sizes of small (0.5 mm in diameter), medium (1 mm in diameter) and large (1.5 mm in diameter). Inspired from animal models, different distances between initial tumor and parent vessels (1, 1.5, 2, 2.5, and 3 mm) are examined. Also two different numbers of initial sprouts (three and five sprouts) are tested. It is assumed that fibronectin is initially secreted by ECs and its concentration is higher near vessels and less near the tumor site. Initial profile of fibronectin concentration is defined as equation (28) 42 : where e 1 and k f are 0.45 and 0.75, respectively. VEGF is distributed within the tissue primarily by quiescent core of the tumor. Since the equations governing angiogenesis are solved in non-dimensional format, values for VEGF and fibronectin concentrations are normalized. Initial VEGF concentration is correlated to initial number of quiescent tumor cells. Initial concentrations of fibronectin and VEGF for different conditions are presented in Supplemental Figures S1 and S2, respectively.

Results
The coupled model of tumor growth and angiogenesis in a vascular tumor microenvironment is established to explore how configuration and properties of parent vessels and initial tumor size and distance regulate tumor growth and angiogenesis. To examine differences among different conditions, the evolution of tumor growth and angiogenesis are presented in 90 days. First, the tumor with following conditions are modeled: initial size of medium (1 mm diameter); the tumor placed at middle of the small domain (4 mm by 4 mm); five sprouts on each parent vessel; and one or two parent vessels placed on walls of the domain (Figure 3(a) and (b)). For the configuration with one single parent vessel on left wall of the domain (Figure 3(a)), branching occurs after 1.3 days and first loop is formed after 1.833 days. The tumor lacks oxygen in first days and all tumor cells turn to hypoxic. Left part of the tumor starts proliferating after 4.5 days and right part of the tumor initiates proliferation after 6.5 days when the capillary network is close enough to the tumor to feed nutrient and oxygen. It should be noted that capillary network reached the tumor in 4 days but capillaries that have reached the tumor site are not functional since they have not anastomosed. Therefore, there is no blood flow in these capillaries. Following propagation of capillary network and formation of closed loops, oxygen diffuses to all parts of tumor and induces proliferation. After a few days of tumor growth when tumor reaches a certain size and consumes oxygen supplied by vessels, farther parts of tumor cannot receive enough oxygen to survive, and therefore they turn to hypoxic cells. The growth of tumor is then biased toward the oxygen source since the right part of the domain does not have access to enough oxygen. Angiogenesis process continues for 7.833 days. The final configuration of tumor growth and angiogenesis after 90 days is shown in Figure 4(a). The results of tumor growth over the course of 90 days for the condition shown in Figure  3(a), with one parent vessel present on left wall and tumor centered in the small domain, are presented in Figure 5 and Supplemental Video 1. The density of tumor cells in initial days of tumor growth increases in the outer rim (where cells are in proliferative phase 2) and cell cycles continue disregarding oxygen availability. As time passes and capillaries provide oxygen for  Through incorporation of one more parent vessel to the right wall (Figure 3(b)), oxygen level in right and left rim of tumor reaches above the hypoxic level and a small portion of tumor starts proliferating at first days. Then, the capillary network propagates and reaches tumor and through preparation of fresh oxygen, all tumor cells start proliferating. Branching initiates after 1 day near right parent vessel. First loop is formed after 2.625 days adjacent to right parent vessel. The other parts of tumor that are in hypoxic state start proliferating after 4 days. The capillary network reaches tumor after 4 days and angiogenesis continues for 7.875 days. The results of 90 days of tumor progression and angiogenesis are illustrated in Figure 4(b). The results of tumor growth and angiogenesis over the course of 90 days are presented in Figure 6 and Supplemental Video 2.
Furthermore, to explore how initial number of sprouts affect tumor growth and angiogenesis, initial number of sprouts is reduced from five to three. Although the simulation process remains the same for each tumor, we study differences in timing of events and explain possible deviations from normal tumor growth discussed above. All information about models as well as references to figures and videos are tabulated in Tables 2 and 3. The conditions for tumor and parent vessels with three sprouts are shown in Figure 3(c) and (d). When there is one parent vessel on the left wall (Figure 3(c)), first branch is observed after 1.125 days and first loop is formed after 2.9 days. However, the tumor remains in hypoxic state until it starts proliferating after 4.5 days when the capillary network supplies enough oxygen. The capillary network reaches tumor after 3.625 days and angiogenesis continues for 6.75 days. As tumor grows, the induced capillary network cannot supply oxygen to right parts of tumor and hence the growth is not eccentric and is directed toward parent vessel and capillaries. Capillary network and tumor growth after 90 days are shown in Figure 4(c). The simulation results of tumor growth and angiogenesis for this configuration over the course of 90 days are illustrated in Supplemental Video 3. When there are two parent vessels with three sprouts on each wall (Figure 3(d)), first branch appears after 0.875 days. First loop is formed after 2.5 days. Other parts of the tumor start proliferating after 3.75 days. Capillaries reach the tumor site after 3.75 days but these capillaries are void of blood flow until they form closed loops. Angiogenesis continues for 6.5 days. Tumor growth and capillary network for this condition after 90 days is shown in Figure 4(d). The simulation results of tumor growth and angiogenesis over the course of 90 days are illustrated in Supplemental Video 4. Comparing the configurations with three and five initial sprouts show that capillary density is reduced for three-sprout configuration, therefore we continue the simulations with parent vessels that have five sprouts.
To understand how tumor size and initial number of cancer cells affect angiogenesis and tumor growth, the model is simulated for a smaller tumor size with the initial diameter of 0.5 mm, with one and two parent vessels and with five initial sprouts. Initial placement of the tumor and parent vessels are shown in Figure 3 Figure  4(f)). The simulation results of tumor growth and angiogenesis over the course of 90 days in this condition are illustrated in Supplemental Video 6.
To investigate how distance between tumor and parent vessel affects tumor growth and angiogenesis, tumors with initial small and medium size are placed closer and farther from the parent vessel which the results are reported in Table 3 for brevity. When tumor with medium size is placed 3 mm away from the parent vessel, the trend of tumor growth is completely different from other configurations (Figure 7(g)). First branches are observed after 5.75 days. First loop is formed after 14.25 days and angiogenesis continues for 26 days. Particularly, right part of the tumor becomes necrotic at first days because of its large distance from oxygen source. Final configurations after 90 days are shown in Figure 8(g), where the results show that tumor growth is directed toward oxygen source in the domain, approaching the newly formed capillary network. Also for the tumor with a distance of 3 mm away from the parent vessel, the duration of angiogenesis increases dramatically and this is due to necrotic region of the tumor with limited production of TAFs (Supplemental Video 13).
The last configuration is to investigate effects of domain size and tumor initial density on vascular tumor growth. First, the centered tumor with an initial medium size and a parent vessel on left wall of the large domain (6 mm 3 6 mm) is simulated (Figure 9(a)). First branch is formed 7 days after initiation of angiogenesis and first loop is formed after 10.5 days. Angiogenesis continues for 30 days. The large size of tumor domain reduces tumor growth due to limited oxygen concentration around the capillary network, and therefore make right part of the tumor necrotic (Figure 10(a) and Supplemental Video 14). Incorporation of another parent vessel on right wall (Figure 9(b)) prevents turning tumor cells to necrotic phase. First branching occurs after 9 days and first loop is formed after 11 days. Second loop is produced after 15 days and it takes 30 days for angiogenesis to form final capillary network (Figure 10 For the large domain (6 mm by 6 mm), the chemotactic strength reduces as the produced VEGF needs to diffuse a longer distance, therefore gradients of VEGF are weak near parent vessels. Since the chemotactic strength within the large domain is limited, the strategy of increasing the tumor density and as a consequence, the number of tumor cells by five folds is applied to make the tumor larger and realize how tumor growth and angiogenesis are impacted. Tumor with initial large size of 1.5 mm in diameter is placed at center of the large domain (Figure 9(c)). For the configuration with a single parent vessel, first branches are observed after 6 days and first loop is formed after 8 days. Angiogenesis also continues for 23 days. Although we enlarged the tumor and rearranged all sprouts to move toward the tumor, the capillary network cannot supply enough oxygen to cells to induce centric tumor growth and fold around the capillary network. The tumor and capillary network are shown after 90 days in Figure  10(c) and illustrated in Supplemental Video 16. Finally, even by adding another parent vessel on right wall (Figure 9(d)), no necrotic part is observed in the tumor site. Branching initiates after 6.5 days adjacent to right parent vessel and first loop is formed after 8.5 days next to right parent vessel. Capillaries reach tumor after 21 days but angiogenesis continues for 28 days. Sprouts on parent vessels could sense VEGF gradient and move toward the tumor. However, the capillary network  could not supply sufficient oxygen to all parts of tumor, therefore a heterogeneous growth of tumor is observed in this configuration (Figure 10(d) and Supplemental Video 17). The number of tumor cells is plotted versus time to compare the trends of tumor growth for different architectures of parent vessels, initial tumor and size of domain (Figure 11(a)). The results confirm that the pattern of tumor growth, when it is fed sufficiently by vessels, obeys the power law of tumor growth, reported by others 55,56 and proved in Ghazani et al. 41 A variety of different growth patterns is observed depending upon their conditions. The diversity of patterns of tumor growth is mainly resulted from gradients of oxygen in microenvironment which dictates biased growth of tumor toward oxygen source and drop of proliferating rim of the tumor. Rapid initial growth rate in first days is due to sufficient space available for cells to grow inside the tumor. When the cells continue proliferation, they are packed inside tumor which increases the density of tumor and enforces the movement toward proliferating rim. Thereby, the rate of tumor growth becomes proportional to the rate of growth in proliferating rim. The larger the tumor diameter, the higher the rate of tumor growth, demonstrated by the increase in slope of growth curve over time (Figure 11(a)). For quiescent state of the tumor, growth is paused for several days but reinitiates once it is fed by new capillaries. This discontinuity is more obvious in asymmetric conditions of the far tumor in the presence of only one parent vessel. In this condition, a necrotic part is formed in the tumor site and the number of tumor cells stays unchanged for several days. However, the growth is recovered after  receiving enough oxygen. Moreover, for the tumor in large domain in the presence of one parent vessel, the slope of growth curve is less than other conditions given the existence of necrotic part and limited oxygen source for tumor cells.
The number of ECs in the domain for 90 days for different configurations of tumor-vessel is shown in Figure 11(b). The trend of variations in number of ECs is identical for all configurations. The rate of increase in number of ECs is slow in first days but increases during sprouting and branching. When the capillary network reaches the tumor, the rate of increase in number of ECs reduces until angiogenesis stops. There is no change in number of ECs after entry of branches into the tumor. Number of ECs is highest in the large domain with two parent vessels. Density of ECs in small tumor configuration is higher than large tumor because ECs have more space near the tumor to branch and anastomose.

Prediction of tumor cells number by artificial neural network
Artificial neural networks have a structure similar to the human brain, which have been widely used in cancer therapy applications. 57,58 The brain, as an information processing system, is made up of key structural elements called neurons. Artificial neural networks also contain a set of interconnected neurons, called layers. 59 To create these layers, these neurons are connected by activation (stimulus) functions. In practice, a limited number of functions are used as activation functions, and neural network researchers often prefer to use nonlinear stimulus functions such as Gaussian or tangent hyperbolic as activation functions. 60 Neural networks, despite their diversity, have a similar structure. A neural network consists of three layers of input, hidden and output layers. The input layer only receives information and acts as an independent variable. Therefore, the Neural networks with radial basis function are widely used to nonparametric estimation of multidimensional functions through a limited set of educational information. Radial neural networks with fast and comprehensive training, are interesting and efficient networks which have received special attention. [61][62][63] Hartman et al. 64 proved that the networks with radial basis function are very powerful approximators. It requires sufficient number of neurons in the hidden layer to make it able to approximate any continuous function with certain accuracy. The interesting aspect is that these networks have this property only with one hidden layer. Radial-based networks derive the most inspiration from the statistical techniques of pattern classification, and their major advantage is the classification of patterns that have nonlinear space. The main RBF architecture consists of a two-layer network such as Figure 12.
The hidden layer establishes a nonlinear adaptation between the inputs and plays an important role in converting nonlinear patterns into linear discriminant patterns. The output layer produces the weighted sum of the linearized patterns along with a linear output. If the RBF is used to approximate the function, such output would be useful; but if patterns need to be classified, then a sigmoid function can be applied to the output nerves to produce output values of 0 or 1. As can be seen from the description above, the unique feature of this network is the process that is performed in the hidden layer. The output of the RBF network can be calculated by equation (29).
Where w j is the weights of each neuron and u j are the centers of the function of each neuron. The famous function in radial networks is a Gaussian or exponential function as follows: In this equation s j is the spread value which is extracted experimentally. The initial size of tumor and distance between tumor and parent vessel are as the inputs of the network and the number of tumor cells is as the output of network. The initial size of tumor, and distance between tumor and parent vessel, calculated in Section 3, have been used as training inputs in the RBF. The schematic of the proposed neural network to predict the tumor volume is represented in the Figure 13. In this figure R is the number of elements in the input vector, p 2 3 1 is the input vector which is made of two elements including initial size of tumor and distance between tumor and parent vessel. Also Q 1 shows  the number of neurons in the hidden layer, Q 2 represents the number of neurons in the output layer, a 1 demonstrates the output of basic radial layer, and a 2 demonstrates the output of RBF networks. Additionally, n 1 is determined by dot product between weights and bias vectors in hidden layer, and n 2 is the summation of weights vectors and bias in output layer. As it is shown in Figure 13, there are nine neurons in the hidden layer and one neuron in output layer. (IW) and (b 1 ) are the values of neuron weights and bias in hidden layer, respectively. Also, (LW) and (b 2 ) are the values of neuron weights and bias in the output layer which are listed in Table 4.
By substituting the above values in equation (31), for every input value (p), the network is able to predict the output values as follows: In this equation, j represents the j-th row of weight matrix in the hidden layer, k shows the k-th element of weight vector and the input vector, i is the number of elements in the weight vector, and n is the number of neurons in layer 1. Also, in this network the value of spread is considered to be s = 1. By considering the above equation and substituting any desired value for tumor distance and initial tumor size in the specified range of acceptable inputs as described in Table 4, the number of tumor cells can be estimated after 90 days of growth. Hence, in order to evaluate the accuracy of the predicted values by equation (31), four arbitrary inputs are given to the network as test inputs and the results of the network are compared to the corresponding results of the mathematical model. The error of the RBF network is demonstrated in Table 5.
The results indicate the high accuracy of the network which is 95%. Using the proposed neural network, the number of tumor cells are predicted at different distances of tumor and the parent vessel in Figure 14. The distance is varied between 1 and 2.5 and the curves for both small and medium tumors are depicted.
According to these results, increasing the distance between tumor and parent vessel leads to a decrease in the number of tumor cells. Also, the rate of decrement for the number of tumor cells in short distances between (1 and 1.5 mm) is less than longer distances (1.5-2.5 mm).

Discussion and conclusion
A new multiscale age-structured cell cycle model coupled with a discrete agent-based angiogenesis method is proposed. Four different phases of tumor cells were incorporated to model cell cycle, including two proliferative phases, quiescent phase and necrotic phase. Different tumor-parent vessel arrangements were simulated in configurations where initial tumor has grown to such a certain size that existing vasculature is unable to provide required oxygen. The scenario demonstrates the situation when tumor induces angiogenesis with the focus on dynamic features of tumor cells and capillary network, population, and density of tumor cells, and sprouting and angiogenesis. In addition to mathematical modeling, an artificial neural network is adopted to be trained for available results to propose an overall formulation for the mathematical model. This formulation would be useful in determining the final number of tumor cells in the domain. As stated in Young, 65 genetic mutations are not the only reason behind tumor cells growth. Many environmental conditions affect tumor growth and can strongly influence the growth of a tumor in oncoming days. With the knowledge of the trend and response of a tumor in subsequent days of growth, it would be possible to suppress and hinder tumor growth. Anti-angiogenesis therapeutics have attracted much attention in recent years. Angiogenesis as the pivotal part of growth of a solid tumor is being studied whether experimentally or numerically. Clinicians are seeking a way to hinder angiogenesis in order to reduce nutrients of the tumor and suppress the growth and metastasis. Prediction of behavior of a tumor in subsequent days in different positions in relation to the tumor will help clinicians better manipulate the therapeutic implementations.
In this study, effects of density of a tumor prior to angiogenesis, number of initial sprouts, and distance of the tumor and parent vessels are investigated. Moreover, to understand the effects of size of the tissue where tumor appears and also in vitro assay cultures, two different domain sizes are studied.
It is demonstrated that the smaller the initial size of the tumor, the stronger the gradient of angiogenic factor is needed to produce sprouts directing straight toward the tumor. The rate of tumor growth for small tumor is less than the rate for large tumor, probably due to a weaker strength of angiogenic chemoattractant gradients. It is also shown that the longer the initial distance of tumor cells to parent vessels, the higher the density of initial tumor is needed to induce sprouting at distant points of vessels.
Even for the configuration of one single parent vessel next to small tumor, if the number of sprouts is adequate and capillary network forms properly, part of the tumor that is not faced to the vessel may receive enough nutrients and oxygen to grow in proliferation state unconditionally. Also the smaller the distance of the tumor to the vessel, the more centric the growth of the tumor.
Approved by many animal models and studied by Folkman, 40 it is stated that when a tumor placed in distances far away from a parent vessel, for example, 2.5 and 3 mm away, it may take much longer time for the vascular network to fully grow and reach the tumor. In this study, it is also concluded that when a tumor is moved in the distances of 2.5 and 3 mm away from a parent vessel, not only a bigger tumor is required to induce sprouting and subsequent angiogenesis phenomenon but also it will take much longer time in comparison to closer distances to tumor to form a vascular network of blood vessels including blood flow.
The RBF neural network provides a comprehensive formulation to be used in prediction of final tumor status in the domain. Based on this formula, the curves for the final number of tumor cells as a response of a tumor in the domain are plotted against tumor-parent vessel distance. As the results demonstrate, beyond a certain distance, both small and medium tumors weaken and the intensity of tumor cell number decrement increases. This gives us the sense that a tumor is strongly dependent upon vascular network such that in farther distances where vasculature growth is delayed, tumor becomes smaller.
All in all, as it is stated in Perfahl et al. 37 and Grogan et al., 38 the geometry of domain in which tumor grows is pivotal and this new hybrid model enables us to monitor growth of a tumor and predict final architecture of the network in different circumstances. Anti-angiogenesis therapeutics are proved to be effective modalities in preventing devastating growth of a tumor and our model again pinpoints the momentous role of angiogenesis in fatal growth of a tumor. Future studies can enrich our model by introducing dynamic tumor network, lymphocytes, stromal cells, and other nutrients such as glucose into the model. Moreover, since blood vessels in tumor-induced neovasculature is tortuous and badly vascularized, smeared representation of neovasculature can be a tool to model patient-specific angiogenesis models which can be combined with the present model to achieve better outcomes for clinical purposes in the future studies.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.

Supplemental material
Supplemental material for this article is available online.