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

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.


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][11][12][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 onedimensional 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.

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,

The Mechanism for the Degradation of Fibronectin
The decay of fibronectin ( ) via protease is assumed to satisfy a reaction mechanism of the form [8], which shows that acts as a catalyst to convert the fibronectin into products .

Modeling Assumptions
In the − plane we envisage a capillary segment of length ℓ 2 microns located along the -axis on the interval [0, ℓ 2 ] with a tumor source located somewhere along the line = ℓ 1 ( Figure 1) [8]. We rescale by /ℓ 1 and by /ℓ 2 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.

The Model
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: where , ] are kinetic parameters, 0 , 0 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: 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 writê =̂( , ) .
A simple transition probability which reflects the influence of enzyme and fibronectin on the motion of endothelial cells iŝ( , ) = 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 takê where Here the , , ( = 1, 2) are positive constants such that 0 < 1 ≪ 1 < 2 and 1 > 1 ≫ 2 > 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 (̂( , ) ) = 0 (at = 0, 1) and take since we assume that the capillary is initially in a rest state.

Solutions of the Linearized Model
We first note that we can write (5) as We look for the linearized solutions of our model using a perturbation of the initial data in the form ( , 0) = 0 = 1, V( , 0) = ( ), ( , 0) = 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] If we now plug (14) into (4) and (11), respectively, we obtain the following equations If one lets → 0 in (15), our linearized model becomes One easily obtains from (16) that If we plug this in (17), we obtain which is a linear equation in the variable . The solution to this equation with the initial condition ( , 0, ) = 0 is If we now write these variables and in (18), it reads For sufficiently large one obtains where = 1 + 2 / . Using the new variables our boundary conditions in (9) become 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 = 0 (at = 0, 1) , and we come up with the following initial-boundary value problem: We now let ( , , ) be any known function and [21] V ( , , ) = ( , , ) − ( , , ) .
If we plug (28) in (27), we obtain We now choose ( , , ) = ( ) to make (29) homogenous. Therefore, since (0) = (1) = 0, (29)-(30) become V ( , 0, ) = − ( ) , 0 < < 1, It is known that this problem has a series solution of the form Letting = 0 in (34) and using (32) give where Therefore, the solution to the problem given by (27) becomes where 0 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: Again, all of the constants and parameters seen in (38) are the same as those mentioned before.

Numerical Example
For numerical purposes we take ( ) = 2 (1 − ) 2 . By the aid of (13) one has = 30. Also we find 0 = −2 and On the other hand, from (20) and (22) 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 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) Also, if we set = 0 in (24), we obtain the steady-state solution as 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.

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.
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 V( , 0) is large and they start with a unimodal distribution in V( , 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 V( , 0) = ( ) = 2 (1 − ) 2 . Even though this function is of unimodal distribution, gradient V( , 0) is not sufficiently large to develop the bimodality of . One could, of course, choose V( , 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 (V, , , ) 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.