Modelling Tumor-induced Angiogenesis: Combination of Stochastic Sprout Spacing and Sprout Progression.

BACKGROUND: Angiogenesis initiated by cancerous cells is the process by which new blood vessels are formed to enhance oxygenation and growth of tumor.
OBJECTIVE: In this paper, we present a new multiscale mathematical model for the formation of a vascular network in tumor angiogenesis process.
METHODS: Our model couples an improved sprout spacing model as a stochastic mathematical model of sprouting along an existing parent blood vessel, with a mathematical model of sprout progression in the extracellular matrix (ECM) in response to some tumor angiogenic factors (TAFs). We perform simulations of the siting of capillary sprouts on an existing blood vessel using finite difference approximation of the dynamic equations of some angiogenesis activators and inhibitors. Angiogenesis activators are chemicals secreted by hypoxic tumor cells for initiating angiogenesis, and inhibitors of the angiogenesis are chemicals that are produced around every new sprout during tumor angiogenesis to inhibit the formation of further sprouts as a feedback of sprouting in angiogenesis. Moreover, for modelling sprout progression in ECM, we use three equations for the motility of endothelial cells at the tip of the activated sprouts, the consumption of TAF and the production and uptake of Fibronectin by endothelial cells.
RESULTS: Coupling these two basic models not only does provide a better time estimation of angiogenesis process, but also it is more compatible with reality.
CONCLUSION: This model can be used to provide basic information for angiogenesis in the related studies. Related simulations can estimate the position and number of sprouts along parent blood vessel during the initial steps of angiogenesis and models the process of sprout progression in ECM until they vascularize a tumor.

Introduction

Angiogenesis is a physiological process involving sprouting of new vessels from a pre-existing vasculature in response to chemical stimuli [1]. Unlike vasculogenesis which depends upon precursor cells [2] or intussusceptive angiogenesis, which splits blood vessels [3], angiogenesis is the process of forming new blood vessels from the pre-existing ones through the migration and proliferation mechanisms [4,5]. Although angiogenesis occurs in normal physiology, at exercise or in the process of wound healing [6], there is a list of pathological conditions involving angiogenesis. Cancer, vascular disease, stroke, neurodegeneration, diabetes, inflammation, asthma, obesity and arthritis are some of these pathologic conditions.

Angiogenesis initiated by cancerous cells is a process by which new blood vessels are formed to enhance oxygenation and growth of tumors [7]. Initially, solid tumors rely on the diffusion from nearby vessels to supply oxygen and nutrients and to remove waste products [8]. As a tumor grows, the demand of oxygen and nutrients increase until the flux of oxygen through the surface of the tumor is too small to supply the entire mass of the cells [8]. A necrotic core of dead cells develops at the center and eventually the tumor stops growing and reaches a steady state size of 2-3 mm in diameter in which the number of dying cells counterbalances the number of proliferating cells [8]. Tumors only grow further if cancerous cells acquire one of the so-called hallmarks of the cancer, the ability to induce angiogenesis through mutations [4,9,10].

Initial Steps of Tumor-induced Angiogenesis

Formation of new blood vessel network is guided by a complex interplay of both pro- and anti-angiogenic molecules produced by a variety of sources including tumor cells, endothelial cells, extracellular matrix, pericytes and plasma clotting products [8,11,12]. The first event of tumor-induced angiogenesis involves the secretion of a number of chemicals, collectively known as tumor angiogenic factors (TAFs) from the cells of a solid tumor into the surrounding tissue [13]. These angiogenesis activators diffuse through tissue space, creating a chemical gradient between the tumor and any existing vasculature [14]. Upon reaching any neighboring blood vessels, a number of chemical interactions between parent blood vessel and pro- and anti- angiogenic factors are done for spacing new sprouts.

Among various anti-angiogenic factors discovered, including tumstatin, arrestin and canstatin, the most notable ones are angiostatin [15-22] and endostatin [18]. It has been shown that both angiostatin and endostatin inhibit sprouting angiogenesis in a dose-dependent manner i.e., the higher the antiangiogenic factor concentration that is present, the more inhibition of angiogenesis is observed [15].

Endostatin is an 18-22 kDa fragment of collagen XVIII, which is notably present amongst other collagens in the vascular and epithelial basement membrane [23]. As tumor cells grow, they degrade the surrounding tissue or ECM in order to spread. Upon this degradation, a number of matrix degradative enzymes (MDEs) such as plasminogen activator (PA) and a large family of matrix metalloproteinases (MMPs) are produced by tumor cells. Another important role of these enzymes is the cleavage of endostatin from collagen XVIII [24]. Karihaloo et al. [25] proposed that MMPs degrade the basement membrane at the tip of the growing ureteric bud and result in an accumulation of endostatin which then acts to inhibit further branching. It has been shown that endostatin also inhibits the endothelial cell proliferation and migration with a possible mechanism to disrupt cell-matrix interactions [12]. Additionally, it has been hypothesized that endostatin plays a role in preventing unchecked outgrowth of the developing ureteric bud and also it acts in a feedback mechanism during angiogenesis [26].

Another inhibitor of the angiogenic sprouting is angiostatin molecule. Angiostatin is a 38 kDa protein, cleaved from serum protein plasminogen by the action of such factors as tissue plasminogen activator (tPA) and several MMPs [27]. It has been shown that during tumor-induced angiogenesis, growth factors secreted by the tumor, initially cause hyper-permeability in the nearby blood vessels [28]. Before spouting is initiated, the breakdown of basement membrane and hyper-permeability in the nearby blood vessels occur and result in the leakage of blood plasma into the surrounding ECM. This plasma contains both plasminogen, from which angiostatin is derived, and the factors which cleave plasminogen to form angiostatin such as tissue plasminogen activator and MMPs [12]. Tissue plasminogen activator and MMPs may also be present in the ECM during angiogenesis through the production by endothelial cells, or may be released by tumor cells. Thus, angiostatin may be formed in the vicinity of blood vessels undergoing angiogenesis [12].

Production of Micro-vessel Structures in ECM

When the interaction between angiogenesis activators and inhibitors activates some endothelial cells along the parent blood vessel, these cells form sprouts. The initial response of sprouts to angiogenic factors is a chemotactic one initiating the migration of cells towards a tumor [29]. The cells continue to make their way through ECM which consists of interstitial tissue, collagen fiber and fibronectin (FN) as well as other components. Among with these components, fibronectin has a significant role in sprout progression. It has been found that endothelial cells synthesize and secrete cellular fibronectin which remains bound to the matrix and does not diffuse [30-34]. The main function of this secreted fibronectin is the adhesion of cells to the matrix and directional movements of a number of cell types [35-37]. Thus, in addition to the chemotactic response of endothelial cells to TAF concentration, there is a complementary haptotactic response to the gradient of adhesiveness of bound fibronectin termed haptotaxis [14]. As endothelial cells migrate towards a tumor, they can form loops and branches and eventually complete angiogenesis, connect with the tumor and penetrate into it to supply the tumor with nutrients that tumor needs to grow further. So, there is also the possibility of tumor cells finding their way into the circulation and being deposited in distant sites in the body resulting in metastasis [38].

In the past several years, a lot of mathematical models have been developed to describe tumor angiogenesis process as an important step in the growth of solid tumors. Some of these models have described the initial steps of tumor angiogenesis including sprout spacing process. Orme and Chaplain [39] developed a three species model of sprout formation that involves endothelial cells, ECM and adhesive sites. They developed their model based on the phenomena occurring within the parent vessel; they used a stability analysis to discuss that a natural spacing appears with the space between sprouts determined from the equations but their model does not include any feedback mechanism [12]. Levine et al. [40] have also developed a model based on the theory of reinforced random walks for the onset of capillary formation. They used the assumption that angiogenic factors are transformed into a proteolytic enzyme that enables a sprout to form. Their model makes predictions about the aggregation of endothelial cells and the perforation of the basement membrane that allows the formation of a new sprout [12]. Again, there is no feedback of the form discussed in the current work. Addison et al. [12] developed a simple mechanistic model with two non-interacting species based on corneal pocket experiments. Their model is able to predict a variety of micro-vessel sprouting patterns representing an inceptive attempt to address the question of how unchecked sprouting of the parent vessel is avoided at the initiation of angiogenesis based on the idea that feedback regulation processes play the dominant role.

Another group of papers has modelled the secondary steps of tumor angiogenesis including models that have focused on sprout progression in ECM and their related interactions with ECM and TAF. One of the models that describe tumor angiogenesis without considering the sprout spacing process, is the model of Anderson and Chaplain [14]. Recently, many continuous-discrete hybrid models [41-49] were developed to investigate tumor angiogenesis and also the blood flow through a network based on the early model of Anderson and Chaplain [14]. In these models, a constant number of active sprouts are initiated from the parent blood vessel at the same time to progress in ECM, something that is different in reality. Also, they have considered the time estimation of tumor angiogenesis process from the beginning of sprout progression instead of the beginning of the activator secretion or the sprout formation.

In this paper, we have developed a simple model to provide a new multiscale mathematical model of the tumor-induced angiogenesis which can estimate the position and number of sprouts along parent blood vessel during the initial steps of angiogenesis to model the process of sprout progression in ECM until they vascularize the tumor. We assumed that initiating of sprouts is based on a hypothesis that sprouting is initiated by a threshold concentration of the activator and anti-angiogenic products synthesized in the immediate vicinity of each new sprout. These products will act locally as a feedback mechanism to prevent the overwhelming formation of new sprouts [12]. Our approach to modelling sprout spacing is different from theirs because our finite difference approximation of the dynamic equations of the activator and inhibitor for finding the position of sprouts along the parent vessel limits the minimum distance between sprouts which allows us to have a more realistic estimation of the time and position of the sprouts. Moreover, we have considered a stochastic radius for the initial inhibitor profile produced around every new sprout. The stochastic radius of the inhibitor not only has more compatibility with the reality but also enables the model to have a better prediction of the sprout numbers. We then use the data obtained from this model including the number, site and time of the activated sprouts to produce micro-vessel structures in ECM. In this paper, we focus on the coupling of these two important models resulting in a better spatiotemporal estimation of tumor angiogenesis. So, just chemical stimuli are investigated and mechanical effects of blood flows in micro-vessel structures are ignored in this paper.

Materials and Methods

Sprout Spacing: Initial Steps of Tumor-induced Angiogenesis

In this model, we focus on two key variables: an angiogenesis activator such as VEGF (vascular endothelial growth factor) which is secreted continuously from the hypoxic tumor cells and an angiogenesis inhibitor such as angiostatin or endostatin which is assumed to be secreted instantaneously at the initiation point of each new sprout. The mathematical model consists of a partial differential equation describing diffusion and decay of the activator in ECM until it reaches the parent blood vessel. There are generally two conditions for every new sprout to form: if the concentration of the activator exceeds a specific trigger level for a point along parent blood vessel, and the concentration of the inhibitor falls below a threshold level for that point, then a new sprout will be formed. As we noted earlier, it has been proposed that MMPs (which are in the vicinity of tumor-induced angiogenesis) degrade the basement membrane at the tip of the growing bud resulting in an accumulation of endostatin which then acts to inhibit further branching [25]. Felbor et al. [26] also suggested that endostatin acts in a feedback mechanism during angiogenesis. Angiostatin may be formed in the vicinity of the blood vessels undergoing angiogenesis too. We assume that, after the formation of each new sprout, depending on the amount of the produced inhibitor in the vicinity of it, the distance between this sprout and the next sprout to form is variable. Since the presence and amount of the inhibitors released by the tumor and endothelial cells are not deterministic, we assume a circular concentration profile with a stochastic radius for the initial inhibitor profile. Each new sprout will increase the concentration of inhibitor in ECM and restricts the positions satisfying sprout formation conditions. For this purpose, we specify the initial inhibitor profile with a circle centered at the site of the new sprout with a random radius of r.

With the exception of the first sprout that is required to satisfy only the condition of the activator for sprouting, all other sprouts need to satisfy both two conditions of the activator and inhibitor concentrations. So, if we denote the angiogenesis activator concentration with A and the angiogenesis inhibitor concentration with I, then, the required conditions for sprouting will be formulated as:

Where, subscripts specify the location on the grid and superscripts specify the time steps. So, if the concentration of the activator around each point of the parent blood vessel becomes greater than or equal to a chosen trigger value, and the inhibitor concentration at that point be less than or equal to a chosen inhibitor threshold value, then, a new sprout will form.

To model activator and inhibitor concentration distributions in the model domain, we have used simple reaction-diffusion equations. We have assumed that the activator is produced by hypoxic tumor cells to simply diffuse. The partial differential equation describing the activator distribution is as follows:

(1)

Where D is the activator diffusion coefficient and λ is its decay rate, both are taken as constants [12].

The inhibitor is produced by every new sprout and then diffuses into the tissue around it to inhibit the formation of further sprouts with a random distance. The equation governing the distribution of the inhibitor concentration is therefore given by:

(2)

Where D is the inhibitor diffusion coefficient and λ is its decay coefficient, both are taken as constants [14].

In order to solve the above equations numerically, the dimensionless parameter values are used for simulations. We have rescaled distance with an appropriate length scale L, time with τ , activator concentration with A and inhibitor concentration with I (where A and I are appropriate reference variables). Therefore, setting

, , , ,

and dropping the tildes for clarity, we have obtained the dimensionless equations,

(3)

(4)

Where:

, , , ,

Then, we have used a finite difference approximation to discretize the above equations as follows:

(5)

(6)

Where, subscripts specify the location on the grid and superscripts specify time steps. l and m are positive parameters which specify the position of variables on the 2D grid i.e. x = lh and y = mh. Time discretization is represented by t = qk where i=1 for activator and 2 for inhibitor equation.

We assumed that the activator and inhibitor remain within the domain under consideration and therefore no-flux boundary conditions of the form:

ξ(-D (7)

for the activator and

ξ(-D (8)

for the inhibitor are imposed on the boundaries of the domain where ξ is an appropriate outward unit normal vector.

As we noted, the first event of the tumor-induced angiogenesis is the secretion of some angiogenesis activators by hypoxic tumor cells. We assumed that hypoxic tumor cells with the radius of 0.1 (in the dimension less area of the simulation) secrete activator. So, if we consider circular geometry for the tumor region, then, the initial condition for the activator has the form,

(9)

In order to find sprout sites, we considered the activator as a continuous source produced by tumor cells which diffuses into ECM and a concentration gradient is established between the tumor and the parent vessel. As we noted before, the newly formed sprout will be an inhibitor source to limit the formation of next sprouts at a stochastic radius of the parent blood vessel.

We assumed initially a circular inhibitor gradient with a random radius of r1 produced around the first sprout. This radius will affect the time and position of the second sprout formation. As soon as the two required conditions 1 and 2 are satisfied for any point of the blood vessel, the second sprout will be formed. This sprout plays the role of an inhibitor source for next sprouts and will increase the concentration of inhibitor in ECM.

As mentioned in [12], we chose the steady state level of the activator at the outermost site that a sprout may form to be the lowest level of activator for which sprouting may occur. This may be used as an activator trigger level. So, the dimensionless expression for the activator steady state may be written in the form of a modified Bessel function of the second kind:

(10)

Values for the inhibitor threshold level may be chosen to give different regimes of the sprouting along a parent blood vessel; a number of them are discussed in the result section of [12]. We have chosen the normalized value of 0.05 for the threshold level of the inhibitor in this paper.

Sprout Progression in ECM

We have used the data obtained from previous stochastic sprout spacing model to determine the site and number of initial activated sprouts, which also considers time differences between sprouts, to initiate sprout progression model.

This model focuses on three key variables: TAF (e.g. VEGF as an important tumor angiogenic factor), fibronectin concentrations and endothelial cells at the tip of the activated sprouts specified in sprout spacing model. This model consists of a system of partial differential equations describing the production and uptake of fibronectin, uptake of TAF by endothelial cells and the chemotaxis and haptotaxis effects of endothelial cells in response to TAF and fibronectin gradients, respectively.

We denoted the endothelial-cell density per unit area by n, TAF concentration by c and fibronectin concentration by f. As already discussed, we assumed that the movement of endothelial cells (at or near a sprout tip) is influenced by three factors: random motility, chemotaxis in response to a gradient of TAF and haptotaxis in response to gradients of fibronectin. Thus, the equation describing changes in endothelial cell density is as follows [14]:

(11)

Where x(c) is a function of the form (x. Also, D and ρ are the random migration, chemotaxis and haptotaxis coefficients, respectively.

As soon as TAF is secreted, diffuses into ECM sets up a concentration gradient between the tumor and any pre-existing vasculature. It should be mentioned that the steady-state of the TAF gradient between the tumor and nearby vessels provides us with the initial conditions for TAF concentration profile. As endothelial cells migrate through ECM in response to this gradient, there is some uptake and binding of TAF by these cells. This process can be easily demonstrated by the following consumption function [14]:

(12)

Where λ is a positive constant to represent TAF consumption rate.

Endothelial cells at the tip of the activated sprouts themselves secrete fibronectin which binds to ECM and does not diffuse. Thus, the equation describing the influence of endothelial cell density on the fibronectin concentration does not contain any diffusion term. Moreover, there are some connections between fibronectin and endothelial cells as they migrate towards the tumor [30]. The production and degradation processes of fibronectin can be modelled by the following equation [14]:

(13)

Where w and μ are positive constants and characterize the production rate of an individual endothelial cell to produce fibronectin and degradation of fibronectin depending on the endothelial cell density, respectively.

Thus, a complete system of equations describing the response of endothelial cells, TAF and fibronectin concentration partially described in the previous parts, reads:

(14)

In order to solve the above system of equations numerically, the dimensionless parameter values are used for simulations. We rescale distance with L (the distance between parent vessel and tumor), time with τ (where D is the TAF diffusion coefficient), endothelial-cell density with n and TAF and fibronectin concentrations with c and f, respectively (where n, c and f0 are appropriate reference variables). So, we obtain non-dimensional system as follows [14]:

(15)

Where D and subject to the no-flux boundary condition,

(16)

on the boundaries of the unit square for endothelial cells. No boundary conditions can be imposed on c and f.

We assumed that tumor is circular and centered at x=1 and y=0.5 with the radius of 0.1 in a dimensionless space (equal to 0.2 mm in the dimensional space). The initial conditions used in the simulations for TAF and fibronectin concentrations are considered as follows respectively: [14]

(17)

(18)

Where v, k and ε are positive constants and r is defined as follows:

(19)

The angiogenesis model incorporates rules for sprout branching and anastomosis and also contains an element of stochasticity for the movement of cells. This discrete model is derived from a discretized form of partial differential equations of the system (15), (using Euler finite-difference methods [50]) and then, the resulting coefficients of the five-point finite-difference stencil are used to generate the movement probabilities of an individual cell in response to its local milieu. The full discretized model is as follows [14]:

(20)

Where, subscripts specify the location on the grid and superscripts specify the time steps. l and m are positive parameters which specify the position of variables on the 2D grid i.e. x = lh and y = mh. Time discretization is represented by t = qk3. The exact forms of P-P involve functions of fibronectin and TAF concentrations near an individual endothelial cell [14]. These coefficients can be thought of as being proportional to the probabilities of endothelial cells being stationary (P) or moving left (P), right (P), up (P) or down (P). The coefficient P, which is proportional to the probability of no movement, has the form:

(21)

and coefficients P, P, P and P, which are proportional to the probabilities of moving left, right, up and down respectively, have the forms:

(22)

(23)

(24)

(25)

Sprout Branching and Anastomosis

Sprout branching as a process by which new sprouts are created from the capillary cells, and anastomosis as a process by which tip cells from one sprout merge with the capillary cells from another sprout are explicitly incorporated to the discrete model. After satisfying three conditions, a capillary sprout can branch at its tip and generate a new sprout. These conditions are as follows [14]:

1. The age of the current sprout is greater than some threshold branching age Ψ (equivalent to a dimensional time of 0.75 days)

2. There is sufficient space locally for a new sprout to form (we assume a forward diagonal branching pattern).

3. The endothelial-cell density is greater than a threshold level n, where n.

If all of the above three conditions are satisfied, we assume that each sprout tip has a probability, P, of generating a new sprout (branching) and this probability is dependent on the local TAF concentration. The sprout branching probabilities associated with various TAF concentration ranges have been chosen on a qualitative basis and are given in Table 1. Figure 1 shows different patterns of sprout branching and anastomosis.

Table 1

Sprout Tip Branching Probabilities as a Function of Local TAF Concentration [14]

TAF Concentration	Sprout Tip Branching Probability
≤0.3	0.0
[0.3-0.5]	0.2
[0.5-0.7]	0.3
[0.7-0.8]	0.4
[0.8-1]	1

TAF: Tumor Angiogenic Factor

Figure1

A Schematic Diagram of Sprout Branching and Anastomosis

As sprouts progress towards a given tumor, driven by the movement probabilities of (20) at each time step of the simulation, endothelial cells at the sprout tips can move to any of the four orthogonal neighbors on the discrete grid. If upon one of these movements, another sprout is encountered, then anastomosis can occur [14].

Results

Simulations of this model were carried out on a 200 × 200 grid, which is a discretization of a unit square [0,1] × [0,1] with a space step of h=0.005 that presents tissue with dimensions [0, 2 mm] × [0, 2 mm]. We assume that tumor (that is the source of angiogenesis activator) is centered at x=1 and y=0.5 with the radius of 0.1 in the dimensionless space (equal to 0.2 mm in the dimensional space) and the parent blood vessel is placed on the left border. Figure 2 shows a schematic diagram of the model domain which includes the tumor on the right boundary of the domain, endothelial cells lining the parent blood vessel on the left boundary and the ECM spaced between them. Discretization of the unit square with a space step of h =0.005 according to our length unit of 2 mm, means h is equivalent to a dimensional length of 10 μm, i.e., approximately the length of one or two endothelial cells [14,51] which has also been schematically shown in Figure 2.

Figure2

A Schematic Diagram of the Model Domain including Parent Blood Vessel, Endothelial Cells Lining Blood Vessel and Tumor as the source of Activator. Discretization of Simulation Domain is also illustrated.

As we discussed before, we choose the time scale for diffusion of the activator and progression of the micro-vessels to be τ and τ, respectively. An estimate for the diffusion coefficient of the activatorD, is 0.864 mm based on the research by Ambrosi et al. [52] using Einstein-Stokes formula [12] and we take D for simulations of sprout progression based on [14]. So, the time scale for the sprout spacing model is equivalent to about 1.15 days and for the sprout progression is 1.59 days. We choose γfor decay coefficient of the activator. The activator decay coefficient plays an important role in determining the patterns of angiogenic sprouting. Estimation for value of λhas been calculated from an in-vitro half-life of VEGF found by Serini et.al. [53] to be 64± 7 min [12]. Choosing the value of 1.15 days for τ1 results the value of 18.08 for γ (details are presented in Table 2) but behavior of the model for a number of decay rates around this value has been investigated and discussed in the following.

Table 2

Summary of Model Parameters

Parameter	Description	Value1	Equation	Reference
β2	Fibronectin Production	0.05	(15)	[14]
D1	Diffusion Coefficient of Activator	0.25	(3)	*
D2	Diffusion Coefficient of Inhibitor	0.25	(4)	**
D3	ECs Random-motility Coefficient	0.00035	(15)	[15]
γ	Fibronectin Uptake	0.1	(15)	[15]
γ1	Activator DecayR	18.08	(3)	***
γ2	Inhibitor Decay Rate	5	(4)	****
ρ	ECs Haptotactic Coefficient	0.34	(15)	[14]
η	TAF Uptake Coefficient	0.1	(15)	[14]
α	Chemotactic Sensitivity Coefficient	0.6	(15)	[14]
*,** Normalization of activator and inhibitor dynamic equations with space and time scales h and τ1, respectively results in values of D1 and D2
D1= (Da τ1) ⁄ L2 = 1/4	D2=(Di τ1) ⁄ L2 =Di ⁄ 4Da =β ⁄ 4 = 1/4
***, *** By definition, γ1 and γ2 are dimension less coefficients of the activator and inhibitor, respectively which are determined based on the type of molecule and its half-life (t1/2):
γ1= τ1 λa=γVEGF=τ1×1 ⁄ (1.44 × t1 ⁄ 2)
γ2= τ1 λi= γ angiostatin=τ1×1 ⁄ (1.44 × t1 ⁄ 2)
In this paper, we consider VEGF molecule as an angiogenesis activator. The half-life of VEGF found by Serini et.al. [53] to be 64± 7 min and the half-life of angiostatin is considered between 4-6 hour based on the paper by Anderson and Chaplain [14].

All values are dimensionless.

All numerical solutions presented in this section were obtained from the finite difference approximations of equations (5) and (6) for sprout spacing model and system (15) with the boundary and initial conditions (16) and (17)-(19) for the sprout progression model, respectively. The iterative steps for the numerical simulations are as follows:

1. Set the boundary and initial conditions for sprout spacing model as given in equations (7)-(9), respectively.

2. Each time step of the simulation process involves solving discrete equation of (5) numerically to obtain the values of the activator at each point of the simulation network and near parent vessel.

3. Check out the first condition of sprouting for each point of the parent vessel to find the site and time of the first sprout formation. If this condition is satisfied, the first sprout will form.

4. Generate r, a random number between 0 and 1, and then initialize the inhibitor instantaneous source with a random radius of r centered at the site of the first sprout.

5. Check out the two conditions of sprouting for either sides of the site of the first sprout. Depending on the radius of r and TAF diffusion and decay coefficients, position and time of the next sprout formation will change.

6. If the two conditions of sprouting satisfy for a point of blood vessel, then, a new sprout will form at that point. After generating a new random number between 0 and 1 as a new radius of r, add the newly produced inhibitor instantaneous source with new radius of r to the existing inhibitor concentration in ECM from previous sprouts.

7. Repeat steps 5 and 6 until there is not any point along the parent blood vessel satisfying sprouting conditions. It is clear that this time is changeable and it is proportional to different values of r produced for the inhibitor sources in previous steps.

8. Set the boundary and initial conditions for sprout progression model as given in equations (16) and (17)-(19) for TAF and fibronectin, respectively. Using the data obtained from sprout spacing model to initialize the number, site and time of initial activated sprouts.

9. Each time step of sprout progression simulation involves solving discrete system of (20) numerically to obtain the values of f and c and then generate five coefficients P–P according to equations (23)-(25). Each endothelial cell is therefore restricted to move to one of its four orthogonal neighboring grid points or remain stationary at each time step.

10. Probability ranges are then computed by summing coefficients P–P to produce five ranges, R=[0, P] and , where j = 1-4.

11. Then, generate a random number between 0 and 1 and, depending on the range into which this number falls, the current individual endothelial cell under consideration will remain stationary (R) or move left (R), right (R), up (R) or down (R).

12. Check out the conditions for sprout branching or anastomosis in each time step.

The parameter values used in the following simulations are dimensionless and given in Table 2. We assumed that the steady state equation of TAF establishes TAF gradient between tumor and the nearby vessels to provide us with the initial condition for TAF concentration profile as given in equation (17). Figure 3 shows the initial conditions for the activator, TAF and fibronectin concentrations in sprout spacing and sprout progression models, respectively.

Figure3

Initial Conditions of (a) the Activator for Sprout Spacing Model, (b) TAF and Fibronectin for Sprout Progression Model.

We start the simulation time from which the hypoxic tumor cells have grown to the radius of 0.1 (dimensionless space) and they have secreted the angiogenesis activator in response to hypoxia. This assumption does not necessarily mean all other tumor cells have accumulated in this radius. So, the total tumor mass may have reached the radius of 1mm in this situation but we just consider TAF secreted by hypoxic tumor cells.

Figure 4 shows the spatiotemporal distribution of the activator in ECM for the value of γ=18.08. As we have illustrated in this figure, the concentration of the activator around a parent vessel is very low in the first steps of simulation and gradually it increases and reaches the activator trigger level at the time of the first sprout formation.

Figure4

Diffusion of Activator Produced by Hypoxic Tumor Cells in ECM to Stimulate the Sprouting of the Blood Vessel on the Left Boundary of the Domain

Figure 5 shows an example of the result of sprout spacing simulation for finding the position of sprouts based on parameter values given in Table 2. In this example, at the end of the simulation, six sprouts will form along parent blood vessel. The first sprout will form about at t =12.81 days and x=0, y=0.51 and it produces a stochastic inhibitor source with the radius of r=0.13. So, the next sprout will distance 0.13 (0.26 mm) from the first sprout in this example. As we observe in Figure 5, the second sprout will form at t =13.15 days and x=0, y=0.38 and produce an inhibitor source of radius r=0.098. Following sprouts will form at different times and produce an instantaneous source of inhibitor to prevent formation of new sprouts around them with a stochastic radius of r.

Figure5

Formation of Sprouts along the Parent Vessel as a result of Numerical Simulation of Sprout Spacing Model with Parameter Values given in Table 2. Six sprouts will form finally in this example.

In the next simulation, we have investigated the effect of activator decay rate, γ, on the behavior of the model with changing γ to be 15 and 20, respectively. Figure 6 shows an example of sprout spacing model considering γ=15. In this case, there are eight sprouts to be activated at the sites of 0.5, 0.96, 0.57, 0.11, 0.37, 0.8, 0.09 and 0.97, respectively. In this regime, it needs more time for the formation of the first sprout than the one we observed in Figure 5; it takes about 18.98 days to form all sprouts.

Formation of Sprouts along Parent Vessel as a result of Numerical Simulation of Sprout Spacing Model with Parameter Values given in Table 2 except γ=15. Eight sprouts will form finally in this example.

Figure 7 shows an example of sprout spacing result considering γ=20. As we observe in this regime, a cluster of sprouts are initially activated in the middle of the domain. All sprouts in a cluster will be formed immediately behind each other at the same time. Formation of these clusters emerging in some activator decay rates is due to our approach of modelling i.e. space discretization. The number of sprouts in these clusters depends on the sharpness of activator profile as well as the space discretization distance. As we have shown in Figure 8, the sharper profiles of the activator or the longer space discretization distances result in fewer numbers of sprouts in a cluster. As we have estimated from the results of simulations, values greater than 20 for the activator decay rate, γ, result in the formation of such clusters. Therefore, changing the type of activator molecule under consideration (VEGF in this paper) will change the value of γ according to its half-life and results in relative changes in sharpness of the activator profile. Choosing the value of space scale greater than h=0.005 not only reduces the precision of modelling, but also decreases the number of sprouts in the cluster.

Figure7

Formation of Sprouts along Parent Vessel as a result of Numerical Simulation of Sprout Spacing Model with Parameter Values given in Table 2 except γ=20. Nine sprouts will form finally in this example.

Figure8

Different Patterns of Sprout Spacing as a result of Numerical Simulation of Sprout Spacing Model for different Values of Activator Decay Rate γ=10, 15, 18.08, 20, 25 and 30, respectively. As we observe, formation of sprout clusters are evident for values of γ≥20.

As we observe in Figure 7, the first three sprouts are formed at the same time in a cluster. We assume that in this situation, each sprout produces a separate resource of inhibitor with a specific radius of r, but because they form at the same time, their inhibitor sources do not affect each other to limit the distance of the next sprout. The fourth sprout forms at time t =13.23 days in the site of 0.23 and the last sprout will form at t =15.12 days in the site of 0.01.

Figure 8 shows various patterns of sprouting with different values of activator decay rates, γ=10, 15, 18.08, 20, 25 and 30, respectively. As we see, the number of sprouts and also the time of the first sprout formation vary with different values of γ.

Figure 9 shows a schematic diagram of different numbers of sprouts specified in a cluster as the first sprouts are formed in the sprout spacing model.

Figure9

A Schematic Diagram of Different Numbers of Sprouts Specified in a Cluster as the First Sprouts are Formed in the Sprout Spacing Model. The sharper profiles of the activator (subplots of (c) and (d)) or the longer space discretization distances (subplots (b) and (d)) result in fewer numbers of sprouts in a cluster.

Sprouts specified in sprout spacing model are used to initiate endothelial cells at the tip of sprouts in the sprout progression model. As we explained in introduction, plasma fibronectin diffuses from the parent vessel in response to the basement membrane degradation of the blood vessel at the time of sprout formation. Numerical simulations of sprout progression model with the given equations for TAF, fibronectin and endothelial cells at the tip of the activated sprout results in Figure 10 to Figure 12 for some specified values of the activator decay rates.

Figure10

Spatiotemporal Evolution of Sprouts from a Numerical Simulation of the Coupled Model with Parameter Values given in Table 2.

Figure11

Spatiotemporal Evolution of Sprouts from a Numerical Simulation of the Coupled Model with Parameter Values given in Table 2 except γ=15.

Figure12

Spatiotemporal Evolution of Sprouts from a Numerical Simulation of the Coupled Model with Parameter Values given in Table 2 except γ=20.

Figure 10 shows the spatiotemporal evolution of micro-vessel structure for the value of γ=18.08. As we explained before, a tumor that is a source of activator is centered at the right border of the domain and endothelial cells lining the parent blood vessel are at the left border. As shown in Figure 5, the first sprout is formed at about t =12.81 days. After 0.34 days, the second sprout will form slowly and they migrate toward the middle of the domain where the source of TAF is situated. The third and fourth sprouts will also form at times t =13.2 and 13.85 days, respectively. We also observe some degree of sprout branching and anastomosis at t =13.85 days. These four sprouts penetrate into the tumor by time t = 14.86 days, so they play a more important role in acquiring oxygen and nutrients for hypoxic tumor cells than two other sprouts to form later. The last two sprouts are also oriented towards the center and micro-vessel network are formed completely by the time t =19.13 days.

Figure 11 shows a numerical simulation of the coupled model which includes sprout spacing and sprout progression models for parameter values given in Table 2 except γ=15. By decreasing the activator decay rate γ, a greater value of activator trigger results, and in similar conditions a greater value of activator trigger level needs more time to satisfy the sprout formation conditions. In this regime, the first sprout will form after about 15 days and the tumor will be vascularized by the time t =17.26 days and the micro-vessel network reaches a stable situation at about t =21.35 days. As expected, the model is capable of reproducing ‘brush border’ effect, the fact that there is an increased frequency of branching at the edge of network as capillary sprouts become closer to the tumor [14].

In Figure 12, we can observe the produced capillary network with the sprouts specified in Figure 7. As we noted earlier, choosing the values of greater than or equal to 20 for the activator decay rate γ, results in the production of some sprout clusters. A cluster of sprouts consists of two or more activated sprouts formed behind each other at the same time. One of the important effects of these sprout clusters is the early creation of loops in the presence of these clusters. These loops occurring in ECM can bypass nutrients and so do not have significant effect on tumor growth. More importantly, these loops can also bypass drugs delivered into the tumor. Therefore, it is one of the important reasons of drug treatment failure. We can see the formation of loops at the first steps of micro-vessel formation in Figure 12.

In the last part, we have investigated the length of micro-vessels produced during angiogenesis process. This is an important parameter affecting the efficiency of drug treatments. The short length of micro-vessels, as they are close to the tumor, results in the high number of them around such a tumor that will dilute drugs delivered to it. At first, we assume all of the vessels initiating at the same time and regret the time differences between sprout formations. In other words, we use the number and site of micro-vessels obtained from sprout spacing model in order to initiate the sprout progression model with the same sprout formation times and then investigate the total trend of the micro-vessel length as they progress to the tumor. Figure 13 shows the length of micro-vessels with respect to the number of them for γ=18.08. As we observe from these plots, the length of micro-vessels decreases almost regularly as the number of them increases. In other words, the newly formed micro-vessels have shorter length than the old ones in general. Figure 14 and 15 show the same plot for γ=18.08 and γ=20, respectively but also when we consider the time differences between initiating sprouts in the sprout progression model. As we see in this regime, there are some peaks that divide the plots into several parts in which the length of micro-vessel decreases as the number of branches increases. There are also some irregularities, too. These irregularities that happen when we consider time differences between sprout formation times, not only make the network more complex and indecisive but also show sprout formation during next times. These long vessels are considered as the main vessels and we need them in order to deliver chemotherapy drugs through for a successful drug treatment. We also observe that the total number of micro-vessels increase significantly when we consider the time differences between sprout formations. Indeed, when sprouts move together through ECM, the probability of anastomosis between them increases and therefore the number of micro-vessel branches decreases. As we see in Figure 13, about 1200 micro-vessels at the end of the simulation time are produced while this number is about 3000 in Figure 14 and 15.

Figure13

Length of Micro Vessels for γ=18.08 when all Sprouts Initiate at the Same Time

Figure14

Length of Micro Vessels for γ=18.08 Considering Time Differences between Sprout Formation

Figure15

Length of Micro Vessels for γ=20 Considering Time Differences between Sprout Formation

Figure 15 shows the histogram of micro-vessel numbers in ECM and different parts of a tumor with respect to their lengths. In this figure, we have assumed that sites of the sprouts obtained from sprout spacing model for γ=18.08 but all these sprouts initiate at the same time and different layers have the radius of 0.1, 0.2 and 0.5 from the center of the tumor. We have considered five groups for sprout length based on the maximum micro-vessel length estimated in simulations and show the center of each group on horizontal axis. For a better comparison, we heave omitted the longest 10% of micro-vessels to show the quality of micro-vessels distribution. There were about 1446 micro-vessel pieces in which 2 vessels had the maximum length of 2490 μm in this simulation. So, less than 0.12% of vessels denoted the most length to themselves. According to Figure 16 and its relative values in Table 3, in the most internal layer of a tumor called “necrotic cell region”, there are just about 60 micro-vessel pieces with short lengths. The number of micro-vessels penetrating in the next layer called “semi-necrotic cell layer” increases and there has found some longer vessels, too. The largest number of micro-vessels can be found in the outermost layer of such a tumor called “well vascular cell region” with about 579 vessel piece having different lengths from the shortest to the longest ones. There are about 100 vessel pieces in the extracellular matrix (ECM) in which their number increases as the capillary length increases. The last row of Table 3 shows the total number of vessels existing in the stated different parts and lengths. We observe that there are a large number of short micro-vessels distributed in the outermost layer of the tumor.

Figure16

Histogram of Micro-vessel Numbers in Different Parts of a Tumor and ECM versus their Length

Table 3

The Number of Micro Vessels in Different Parts of Tumor and ECM

Region	[0-49.8]	[49.8-99.6]	[99.6- 149.4]	[149.4-199.2]	[199.2-249]	Sum
Necrotic Cell Region	59	1	0	0	0	60
Semi-necrotic Cells Region	454	103	5	0	0	562
Well Vascular Cell Region	191	224	102	33	29	579
ECM	0	3	13	28	56	100
Total	704	331	120	61	85	1301

Discussion

In this study, we have investigated the process of tumor angiogenesis divided into two important parts of the sprout spacing and sprout progression. This mathematical model enables us to have a better estimation of spatiotemporal evolution of capillary structures.

We have developed sprout spacing model based on the known model by Addison et al. [12] but we have also modified the initial constant radius profile of the inhibitor produced around every new sprout with a random radius profile due to the stochastic nature of inhibitor production and its concentration in the vicinity of a blood vessel during tumor angiogenesis. One of the interesting advantages of this modification is that the number of activated sprouts along a parent vessel is not a fixed number. Our finite difference approximation of the model is the other important modification enabling us to model the parent blood vessel as a discrete entity which limits the distance between the sprouts with a minimum length of 10μm which is equal to the length of one or two endothelial cells. We removed the symmetry sprouting condition of the blood vessel and let the sprouts to form without any special patterns. Investigation of the micro-vessel structures activating sprouts and their branches is also considered in this model. One of the famous models for sprout progression is the model developed by Anderson and Chaplain [14] that uses a finite difference approach to model micro-vessel structures in tumor angiogenesis. They have assumed a constant number of activated sprouts distributed with a specified distance along the parent vessel. Their model does not include the sprout spacing process. Additionally, our mathematical model describes diffusion and decay of the activator (which is secreting continuously from the hypoxic tumor cells) in ECM until it reaches the parent blood vessel. This consideration enables us to estimate the tumor angiogenesis process from the beginning of activator secretion from hypoxic tumor cells instead of the beginning of sprout progression with fixed number of sprouts.

We also modified orthogonal sprout branching patterns with forward diagonal branching patterns in the sprout progression model to comply with the reality. With this modification, we observed faster progression of sprouts towards the tumor and more realistic structures.

Simulation results of this model demonstrate the ability of this model to produce different patterns of micro-vessel structures during tumor-induced angiogenesis. As we observed, the formation site of the first sprout is approximately fixed while the time of the first sprout formation changes regularly with the chosen value for the activator decay rate γ. On the other hand, sites and times of next sprout formation are not constant and depend on the radius of the inhibitor sources and activator diffusion coefficient. However, it takes about 10-15 days (after secretion of activator from hypoxic tumor cells) for the first sprout to form. Theoretically, finite difference approximation of the model limits the maximum number of sprouts up to 200 along a 2 mm length of a blood vessel under ideal circumstances but according to the results of our simulations, the number of sprouts rarely exceeds the range of 5-10 under real conditions.

Moreover, we tried to estimate the relative number and length of micro-vessels in different parts of a tumor and extracellular matrix. Although we have not included any specific features for these different layers, we can see that there are two different regions in a tumor including the most and the least numbers of micro-vessels. Considering different features of these layers may regulate the estimate values for number and length of micro-vessels.

This model can be used as a supplement model to provide mathematical analysis in future studies. In this paper, we focused on the initial steps of tumor angiogenesis. This model can be developed to model blood flow into capillary structures and to investigate angiogenesis under the effects of some mechanical stimuli like shear stress. This model is able to predict the number of micro-vessels in a more realistic situation and also estimate the main blood vessels in produced capillary networks.