Solutions of a Linearized Mathematical Model for Capillary Formation in Tumor Angiogenesis: An Initial Data Perturbation Approximation

Abstract

We present a mathematical model for capillary formation in tumor angiogenesis and solve it by linearizing it using an initial data perturbation method. This method is highly effective to obtain solutions of nonlinear coupled differential equations. We also provide a specific example resulting, that even a few terms of the obtained series solutions are enough to have an idea for the endothelial cell movement in a capillary. MATLAB-generated figures are provided, and the stability criteria are determined for the steady-state solution of the cell equation.

1. Introduction

Angiogenesis is the main feature of neovascularization, the formation of new blood vessels. It is defined as the outgrowth of new vessels from a preexisting vascular network and is fundamental to the formation of blood vessels during placental growth and wound healing, for example. It is also known that it occurs in three sequential steps [1]. First, the endothelial cells (EC) lining the vascular basal lamina (BL) (or basement membrane) degrade this membrane. Second, the EC migrate and proliferate (via mitosis) into the extracellular matrix (ECM). Finally, capillary loops form. In recent years, progress has been made to understand this phenomenon at the molecular level. This includes the identification of potent angiogenic factors, the discovery of the role of proteases, the importance of the ECM, and the emerging characterization of signal transduction pathways in EC.

One of the major components of the ECM is fibronectin, a large, highly adhesive glycoprotein particularly abundant in plasma, connective tissue matrices, and BL [2]. It is also known to enhance EC adhesion to collagen and is produced by EC [3]. The simplest unifying interpretation of these findings is that it functions as adhesive protein that binds cells to other cells or to substrate. Fibronectin-treated cells also migrate more rapidly, both as single cells or as masses of cells migrating out from cell aggregates [3].

As stated in [4, 5] EC are to be stimulated by a tumor angiogenic factor for angiogenesis to occur. Also, active enzyme stimulates the EC migration [6]. Once the EC are stimulated, the long-time tendency of them is towards the transition probability density function (TPDF) [7] of active enzyme and fibronectin (see [8] for mathematical proof of this).

Endothelial cell migration and proliferation also occur during endothelial repair in situ; the ability to penetrate the vascular basement membrane, on the other hand, is an aspect of endothelial cell behavior uniquely expressed during angiogenesis [9].

There have been many mathematical models describing tumor angiogenesis (see [10–13] and references therein). For example, in [11] the authors propose a review and critical analysis of the asymptotic limit methods focused on the derivation of macroscopic equations for a class of equations modeling complex multicellular systems by methods of the kinetic theory for active particles, and in [12] the authors deal with the derivation of macroscopic equations of biological tissues for a class of nonlinear equations, with quadratic type nonlinearity, modeling complex multicellular systems. Also, in [13] a continuous model for three early stage events in angiogenesis, initiation, sprout extension, and vessel maturation, is presented, and in [8, 14], a mathematical model for capillary network formation is presented, and a mathematical analysis of it is given in one and two dimensions, respectively. The mathematical analyses there were very useful to understand mathematically the mechanism of the cell movement in a tissue (see the biological overview in there). On the other hand, in [15] the authors study some qualitative properties of the solutions of a nonlinear flux-limited equation arising in the transport of morphogens in biological systems, and in [16], to eliminate nonbiological behaviors from diffusion models the authors introduce flux-limited spreading, which implies a restricted velocity for morphogen propagation and a nonlinear mechanism of transport.

The layout of the rest of the paper is as follows. First, we describe our modeling assumptions, and write down our one-dimensional model equations originally presented in [17] in detail. Second, we linearize our model using an initial data perturbation method and solve it. Third, we provide a specific example and plot the results. Finally, we close the paper by presenting the conclusions and discussions on our analysis and the biological importance of our results.

2. The Mechanism for the Production of Protease

If 𝑉 denotes a molecule of angiogenic factor (substrate) and 𝑅_𝐸 denotes a receptor on the endothelial cell wall, they combine to produce an intermediate complex, 𝑅_𝐸𝑉 which is an activated state of the receptor that results in the production and secretion of proteolytic enzyme, 𝐶, and a modified intermediate receptor 𝑅^′_𝐸. The receptor 𝑅^′_𝐸 is subsequently removed from the surface by a mechanism that is presumed to be very fast in the time scale of the production of protease 𝐶. The receptor 𝑅^′_𝐸 is then either recycled back to the cell surface to again become 𝑅_𝐸 or degraded and new 𝑅_𝐸 is synthesized, which then moves to the cell surface to replace the 𝑅_𝐸 that had been removed [8].

The point of view is that the receptors at the surface of the cell function in the same way an enzyme functions in classical enzymatic catalysis [8]. In symbols, 𝑉+𝑅_𝐸⇌𝑅_𝐸𝑉𝑅_𝐸𝑉⟶𝐶+𝑅_𝐸. (1)

3. The Mechanism for the Degradation of Fibronectin

The decay of fibronectin (𝐹) via protease is assumed to satisfy a reaction mechanism of the form [8], [][]𝐶+𝐹⇌𝐶𝐹𝐶𝐹⟶𝐶+𝐹^′ (2) which shows that 𝐶 acts as a catalyst to convert the fibronectin into products 𝐹^′.

4. Modeling Assumptions

In the 𝑥−𝑦 plane we envisage a capillary segment of length ℓ₂ microns located along the 𝑦-axis on the interval [0,ℓ₂] with a tumor source located somewhere along the line 𝑥=ℓ₁ (Figure 1) [8]. We rescale 𝑥 by 𝑥/ℓ₁ and 𝑦 by 𝑦/ℓ₂ so that this rectangle becomes a unit square. Therefore, we now have 0≤𝑥, 𝑦≤1. Basically, the problem consists of two parts: (i) the dynamics on the 𝑦 axis, namely, in the capillary (1D problem); (ii) the dynamics in the unit square, namely, in the ECM (2D problem). We couple those two dynamics via some boundary conditions (see [17] for details). The 1D problem is the focus of this paper.

Figure 1 

Extracellular matrix (ECM).

5. The Model

We use the following notation for the concentrations of the various chemical species along the capillary wall in 𝜇M (micromoles per cubic liter): 𝑣𝑐𝑓𝜂𝑦,𝑡=angiogenicfactordensity,𝑦,𝑡=proteolyticenzymedensity,𝑦,𝑡=fibronectindensity,𝑦,𝑡=endothelialcelldensity. (3) If we apply the law of mass action to (1)-(2) (see [17] for details), which asserts that the rate of reaction is proportional to the product of the concentrations of the reactants, we obtain the following system of equations: 𝜕𝑣𝜕𝑡=−𝜆𝑣𝜂1+𝜈𝑣𝜂₀,𝜕𝑐=𝜕𝑡𝜆𝑣𝜂1+𝜈𝑣𝜂₀,𝜕𝑓𝑓𝜕𝑡=𝛽𝑓1−𝑓₀𝜂𝜂₀−𝜇𝑐𝑓, (4) where 𝜆, 𝜈 are kinetic parameters, 𝜂₀, 𝑓₀ are some reference numbers for EC and fibronectin, respectively, and 𝛽, 𝜇 are some positive constants.

Furthermore, if we use the theory of reinforced random walk derived by Davis [18], we obtain the EC movement equation (see [17] for the derivation of the following equation) as follows: 𝜕𝜂𝜕𝑡=𝐷_𝜂𝜕𝜂𝜕𝜕𝑦𝜂𝜕𝑦ln𝜏, (5) where 𝐷_𝜂 is a positive constant, the EC diffusion coefficient in the capillary, and 𝜏 is the so-called transition probability function. This function has the effect of biasing the random walk of endothelial cells. In this case, we know that the walk is influenced by the proteolytic enzyme it produces in response to the angiogenic factor that has made its way to the cell receptors, and by the fibronectin in the BL, thus, we write 𝜏=𝜏𝑐,𝑓. (6)

A simple transition probability which reflects the influence of enzyme and fibronectin on the motion of endothelial cells is 𝜏(𝑐,𝑓)=𝑐^𝛾1𝑓^−𝛾2 for positive constants 𝛾_𝑖(𝑖=1,2).

The biological interpretation of this choice is that endothelial cells prefer to move into the regions where 𝑐 is large or where 𝑓 is small.

In order to avoid singularities in ln𝜏 and its derivatives in (5), it is useful to take 𝜏𝑐,𝑓=𝜏₁(𝑐)𝜏₂𝑓, (7) where 𝜏₁(𝑐)=𝛼₁+𝑐𝛼₂+𝑐^𝛾1,𝜏₂𝑓=𝛽₁+𝑓𝛽₂+𝑓^𝛾2. (8) Here the 𝛼_𝑖,𝛽_𝑖,𝛾_𝑖(𝑖=1,2) are positive constants such that 0<𝛼₁≪1<𝛼₂ and 𝛽₁>1≫𝛽₂>0. Clearly, 𝜏 is not singular for small or large values of 𝑐,𝑓 and will approximate 𝑐^𝛾1𝑓^−𝛾2 over a considerable range of these variables [17].

We impose zero flux boundary conditions for the cells in (5) 𝐷_𝜂𝜂𝜕𝜂𝜕𝑦ln𝜏𝑐,𝑓=0at𝑦=0,1 (9) and take 𝜂𝑦,0=𝜂₀>0, (10) since we assume that the capillary is initially in a rest state.

6. Solutions of the Linearized Model

We first note that we can write (5) as 𝜕𝜂𝜕𝑡=𝐷_𝜂𝜕𝜂𝜕𝑦^𝑦𝜏−𝜂^′1𝜏₁𝑐_𝑦𝜏−𝜂^′2𝜏₂𝑓_𝑦, (11) where 𝜏^′₁𝜏₁=𝛾₁𝛼₂−𝛼₁𝛼₁+𝑐𝛼₂,𝜏+𝑐_′2𝜏₂=𝛾₂𝛽₂−𝛽₁𝛽₁+𝑓𝛽₂.+𝑓 (12) We look for the linearized solutions of our model using a perturbation of the initial data in the form 𝜂(𝑦,0)=𝜂₀=1, 𝑣(𝑦,0)=𝜖𝜃(𝑦), 𝑐(𝑦,0)=𝑐₀=0, 𝑓(𝑦,0)=𝑓₀=1. Here 𝜖 is some positive parameter and 𝜖𝜃(𝑦) is some initial disturbance in growth factor concentration normalized in such a way that [19, 20] ¹₀𝜃𝑦𝑑𝑦=1. (13) Correspondingly let us write: 𝑣,𝑐,𝑓,𝜂.𝑦,𝑡=𝜖𝑤𝑦,𝑡,𝜖𝑦,𝑡=𝜖𝜉𝑦,𝑡,𝜖𝑦,𝑡=1−𝜖𝜑𝑦,𝑡,𝜖𝑦,𝑡=1+𝜖𝑛𝑦,𝑡,𝜖 (14)

If we now plug (14) into (4) and (11), respectively, we obtain the following equations 𝜖𝑤_𝑡()=−𝜆𝜖𝑤1+𝜖𝑛,1+𝜈𝜖𝑤𝜖𝜉_𝑡=𝜆𝜖𝑤(1+𝜖𝑛),1+𝜈𝜖𝑤−𝜖𝜑_𝑡=𝛽𝜖𝜑1−𝜖𝜑(1+𝜖𝑛)−𝜇1−𝜖𝜑𝜖𝜉,𝜖𝑛_𝑡=𝐷_𝜂𝜕𝜕𝑦𝜖𝑛_𝑦−(1+𝜖𝑛)𝛾₁𝛼₂−𝛼₁𝜖𝜉_𝑦𝛼₁+𝜖𝜉𝛼₂−()𝛾+𝜖𝜉1+𝜖𝑛₂𝛽₂−𝛽₁𝜖𝜑_𝑦𝛽₁+1−𝜖𝜑𝛽₂.+1−𝜖𝜑 (15) If one lets 𝜀→0 in (15), our linearized model becomes 𝑤_𝑡𝜉=−𝜆𝑤,_𝑡𝜑=𝜆𝑤,_𝑡𝑛=−𝛽𝜑+𝜇𝜉,_𝑡=𝐷_𝜂𝜕𝑛𝜕𝑦^𝑦−𝜌₁𝜉_𝑦−𝜌₂𝜑_𝑦, (16)(17)(18) where 𝜌₁=𝛾₁𝛼₂−𝛼₁𝛼₁𝛼₂,𝜌₂=𝛾₂𝛽₂−𝛽₁𝛽₁+1𝛽₂+1. (19) One easily obtains from (16) that 𝑤𝑦(),𝜉𝑦().𝑦,𝑡,𝜖=𝜃exp−𝜆𝑡𝑦,𝑡,𝜖=𝜃1−exp−𝜆𝑡 (20) If we plug this 𝜉 in (17), we obtain 𝜑_𝑡+𝛽𝜑=𝜇𝜃𝑦1−exp(−𝜆𝑡), (21) which is a linear equation in the variable 𝜑. The solution to this equation with the initial condition 𝜑(𝑦,0,𝜖)=0 is 𝜑𝑦1𝑦,𝑡,𝜖=𝜇𝜃𝛽−exp(−𝜆𝑡)+𝛽−𝜆𝜆exp−𝛽𝑡𝛽𝛽−𝜆. (22) If we now write these variables 𝜉 and 𝜑 in (18), it reads 𝑛_𝑡=𝐷_𝜂𝑛_𝑦𝑦−𝜌₁1−exp(−𝜆𝑡)+𝜌₂𝜇1𝛽−exp(−𝜆𝑡)+𝛽−𝜆𝜆exp−𝛽𝑡𝛽𝛽−𝜆×𝜃^′′𝑦. (23) For sufficiently large 𝑡 one obtains 𝑛_𝑡=𝐷_𝜂𝑛_𝑦𝑦−𝐾𝜃^′′𝑦, (24) where 𝐾=𝜌₁+𝜌₂𝜇/𝛽.

Using the new variables our boundary conditions in (9) become 𝑛_𝑦=𝜌₁𝜉_𝑦+𝜌₂𝜑_𝑦at𝑦=0,1. (25) Since we want the function 𝜃(𝑦) to be of unimodal distribution type, it is reasonable to choose 𝜃^′(0)=𝜃^′(1)=0. Therefore, the conditions in (25) now read 𝑛_𝑦=0at𝑦=0,1, (26) and we come up with the following initial-boundary value problem: 𝑛_𝑡−𝐷_𝜂𝑛_𝑦𝑦=−𝐷_𝜂𝐾𝜃^′′𝑦,𝑛𝑛𝑦,0,𝜖=0,0<𝑦<1,_𝑦(0,𝑡,𝜖)=𝑛_𝑦(1,𝑡,𝜖)=0,𝑡>0. (27) We now let 𝑀(𝑦,𝑡,𝜖) be any known function and [21] 𝑣𝑦,𝑡,𝜖=𝑛𝑦,𝑡,𝜖−𝑀𝑦,𝑡,𝜖. (28) If we plug (28) in (27), we obtain 𝑣_𝑡−𝐷_𝜂𝑣_𝑦𝑦=−𝐷_𝜂𝐾𝜃^′′𝑦−𝑀_𝑡−𝐷_𝜂𝑀_𝑦𝑦,𝑣𝑣𝑦,0,𝜖=−𝑀𝑦,0,𝜖,0<𝑦<1,_𝑦(0,𝑡,𝜖)=−𝑀_𝑦()𝑣0,𝑡,𝜖,𝑡>0,_𝑦(1,𝑡,𝜖)=−𝑀_𝑦(1,𝑡,𝜖),𝑡>0. (29)(30) We now choose 𝑀(𝑦,𝑡,𝜖)=𝐾𝜃(𝑦) to make (29) homogenous. Therefore, since 𝜃^′(0)=𝜃^′(1)=0, (29)-(30) become 𝑣_𝑡−𝐷_𝜂𝑣_𝑦𝑦𝑣𝑦𝑣=0,𝑦,0,𝜖=−𝐾𝜃,0<𝑦<1,_𝑦(0,𝑡,𝜖)=𝑣_𝑦(1,𝑡,𝜖)=0,𝑡>0. (31)(32)(33) It is known that this problem has a series solution of the form 𝑣=𝐶𝑦,𝑡,𝜖⁰2+^∞_𝑛=1𝐶_𝑛𝑒^{−𝑛2𝜋2𝐷𝜂𝑡}cos𝑛𝜋𝑦,0≤𝑦≤1,𝑡>0. (34) Letting 𝑡=0 in (34) and using (32) give 𝑣𝑦=𝐶𝑦,0,𝜖=−𝐾𝜃⁰2+^∞_𝑛=1𝐶_𝑛cos𝑛𝜋𝑦,0≤𝑦≤1, (35) where 𝐶₀=−2𝐾¹₀𝜃𝑦𝐶𝑑𝑦=−2𝐾,_𝑛=−2𝐾¹₀𝜃𝑦cos𝑛𝜋𝑦𝑑𝑦,𝑛=1,2,…. (36) Therefore, the solution to the problem given by (27) becomes 𝑛=𝐶𝑦,𝑡,𝜖⁰2+^∞_𝑛=1𝐶_𝑛𝑒^{−𝑛2𝜋2𝐷𝜂𝑡}𝑦cos𝑛𝜋𝑦+𝐾𝜃,0≤𝑦≤1,𝑡>0, (37) where 𝐶₀ and 𝐶_𝑛’s are the same as those mentioned before.

Consequently, from (14) the perturbation solutions of the system consisting of (4)-(5) together with the initial and boundary conditions are obtained as follows: 𝑣𝑦()𝑐𝑦()𝑓𝑦1𝑦,𝑡=𝜖𝜃exp−𝜆𝑡,0<𝑦<1,𝑡>0,𝑦,𝑡=𝜖𝜃1−exp−𝜆𝑡,0<𝑦<1,𝑡>0,𝑦,𝑡=1−𝜖𝜇𝜃𝛽−exp(−𝜆𝑡)+𝛽−𝜆𝜆exp−𝛽𝑡𝛽,𝜂𝐶𝛽−𝜆0<𝑦<1,𝑡>0,𝑦,𝑡=1+𝜖⁰2+𝜖^∞_𝑛=1𝐶_𝑛𝑒^{−𝑛2𝜋2𝐷𝜂𝑡}𝑦cos𝑛𝜋𝑦+𝜖𝐾𝜃,0≤𝑦≤1,𝑡>0. (38) Again, all of the constants and parameters seen in (38) are the same as those mentioned before.

7. Numerical Example

For numerical purposes we take 𝜃(𝑦)=𝐴𝑦²(1−𝑦)². By the aid of (13) one has 𝐴=30. Also we find 𝐶₀=−2𝐾 and 𝐶_𝑛=0when𝑛=1,3,5,…,1440𝐾𝑛⁴𝜋⁴when𝑛=2,4,6…. (39)

Therefore, we have 𝑣𝑦,𝑡,𝜖=−𝐾1−90𝜋^4∞_𝑛=1cos2𝑛𝜋𝑦𝑛⁴𝑒^{−4𝑛2𝜋2𝐷𝜂𝑡},𝑛0<𝑦<1,𝑡>0,𝑦,𝑡,𝜖=−𝐾1−90𝜋^4∞_𝑛=1cos2𝑛𝜋𝑦𝑛⁴𝑒^{−4𝑛2𝜋2𝐷𝜂𝑡}+30𝐾𝑦²1−𝑦²,0<𝑦<1,𝑡>0. (40)(41) On the other hand, from (20) and (22) one has 𝑤𝑦,𝑡,𝜖=30𝑦²1−𝑦²()𝜉exp−𝜆𝑡,0<𝑦<1,𝑡>0,𝑦,𝑡,𝜖=30𝑦²1−𝑦²(,𝜑1−exp−𝜆𝑡)0<𝑦<1,𝑡>0,𝑦,𝑡,𝜖=30𝜇𝑦²1−𝑦²×1𝛽−exp(−𝜆𝑡)+𝛽−𝜆𝜆exp−𝛽𝑡𝛽,𝛽−𝜆0<𝑦<1,𝑡>0. (42)

As a result, the perturbation solutions of the system consisting of (4)-(5) for this special choice of 𝜃(𝑦) are obtained by plugging (41)-(42) into (14), and they are 𝑣𝑦,𝑡=30𝜖𝑦²1−𝑦²()𝑐exp−𝜆𝑡,0<𝑦<1,𝑡>0,𝑦,𝑡=30𝜖𝑦²1−𝑦²(,𝑓1−exp−𝜆𝑡)0<𝑦<1,𝑡>0,𝑦,𝑡=1−30𝜖𝜇𝑦²1−𝑦²×1𝛽−exp(−𝜆𝑡)+𝛽−𝜆𝜆exp−𝛽𝑡𝛽,𝜂𝛽−𝜆0<𝑦<1,𝑡>0,𝑦,𝑡=1−𝜖𝐾1−90𝜋^4∞_𝑛=1cos2𝑛𝜋𝑦𝑛⁴𝑒^{−4𝑛2𝜋2𝐷𝜂𝑡}+30𝜖𝐾𝑦²1−𝑦². (43)(44) In this example, we take 𝜖=0.05, 𝜆=20, 𝜇=5, 𝛽=10, 𝐾=0.1, and 𝐷_𝜂=0.25 for our computations. In the following the reader finds the pictures of angiogenic factor, proteolytic enzyme, fibronectin, and endothelial cell densities, respectively, obtained from (43)-(44). The pictures are plotted until 𝑡=0.1.

On the other hand, one easily gets from (41) that lim_𝑡→∞𝑛𝑦,𝑡,𝜖=−𝐾+30𝐾𝑦²1−𝑦²=−𝐾+𝐾𝜃𝑦. (45) Also, if we set 𝑛_𝑡=0 in (24), we obtain the steady-state solution as 𝑛𝑦,𝜖=𝐵+𝐾𝜃𝑦, (46) where 𝐵 is an arbitrary constant. Comparison of (45) and (46) gives that this steady state is stable if 𝐵 is taken to be the constant −𝐾, which is −0.1 in this case.

8. Conclusion and Discussion

In this paper we have presented a mathematical model for capillary formation in tumor angiogenesis and solved it by linearizing it using an initial data perturbation method. Even though our results here are almost the same as those obtained in [17, 19], where this model has been solved by a classical explicit method and by the method of lines, respectively, we believe the linearization method is much more easier and effective.

This model is based on the cell biologist’s observations concerning cell movement. In the case of endothelial cells, these are more likely to move places where fibronectin density is low to follow a chemical trail consisting of growth factor and to respond to growth factor by moving into new space created by an enzyme they produce that in turn destroys fibronectin as well as other ECM components [19].

Figures 3, 4, and 5 show the angiogenic factor, proteolytic enzyme, and fibronectin densities, respectively. We have used only four-term expansion of the series in (44) to draw Figure 6. Even with four-term expansion we have obtained what we expect to see for the endothelial cell movement in a capillary. This also shows the effectiveness of our method. Of course, a more stable solution for the cell equation can be obtained by expanding more terms of the series solution.

Figure 1 and Figures 2–6 have been created using WinEdt/MiKTeX and MATLAB, respectively.

Figure 2 

Profile of possible unimodal distribution of 𝜃(𝑦).

Figure 3 

Angiogenic factor density 𝑣(𝑦,𝑡).

Figure 4 

Proteolytic enzyme density 𝑐(𝑦,𝑡).

Figure 5 

Fibronectin density 𝑓(𝑦,𝑡).

Figure 6 

Endothelial cell density 𝜂(𝑦,𝑡).

When running numerical computations, one might face the bimodality properties of endothelial cells as time increases. But, for what conditions this situation holds? This question has been addressed in [19], and there the authors have observed that if gradient 𝜕_𝑦𝑣(𝑦,0) is large and they start with a unimodal distribution in 𝑣(𝑦,0), a bimodal distribution in ln𝜏(𝑐,𝑓) may develop. The bimodal structure of ln𝜏(𝑐,𝑓) in turn leads to the bimodality of 𝜂, the EC distribution, and provides evidence that our model is in agreement with biology since the bimodal EC distribution correspond to the beginnings of the EC lining of the nascent capillary. As the reader remembers, in our numerical example we take 𝑣(𝑦,0)=𝜖𝜃(𝑦)=𝜖𝐴𝑦²(1−𝑦)². Even though this function is of unimodal distribution, gradient 𝜕_𝑦𝑣(𝑦,0) is not sufficiently large to develop the bimodality of 𝜂. One could, of course, choose 𝑣(𝑦,0)=𝜖𝐴𝑦^𝑚(1−𝑦)^𝑚 for some 𝑚>0 and try to find the desired constant 𝑚 which creates the bimodality of 𝜂. According to our numerical experience, finding such constant 𝑚 may not be easy.

We must also mention the importance of the function 𝜃(𝑦) appearing in the initial data. It is a function that initiates the dynamics in the model equations. If it is zero, every variable in the model stays dormant, and no capillary forms for angiogenesis to occur.

On the other hand, we have determined the stability criteria for the steady-state solution of the endothelial cell equation. Although we do not present the pictures in the paper, all of the four variables (𝑣,𝑐,𝑓,𝜂) come to steady state, about 𝑡=25.

We are in the process of expanding the idea used in this paper to our 2D problem (in ECM, see Figure 1). The results obtained here look promising that we will have similar results for cells in ECM, where a capillary branches from a given capillary and lining itself with endothelial cells as it progresses.