Numerical Simultation of Capillary Formation during the onset of Tumor Angiogenesis

Overview

This web page contains reference material for the paper


Numerical simulation of capillary formation during the onset of tumor angiogenesis,


co-authored by M.W. Smiley , H.A. Levine and M. Nilsen-Hamilton. The reference material consists of several graphs that were produced using the algorithm described in this paper. The graphs are snapshots in time of the density and concentration functions that are used in the modeling this biological process. To put these functions and their graphs in a proper context the introduction to the paper is given below. A complete copy of the paper can be downloaded by using the appropriate link below. For convenience in accessing this link and various portions of the web page the following quick links may be used.

Introduction from the paper

According to a hypothesis proposed by J. Folkman (see for example [1]), a tumor can evolve from an avascular state (without its own blood supply) to a vascular state by releasing proteins into the surrounding medium that stimulate the outgrowth of new capillaries from existing capillaries. Vascularization occurs when the newly formed capillaries reach the tumor. The proteins emitted by the tumor are known as growth factors, and the outgrowth of new capillaries from existing capillaries is known as angiogenesis. A mathematical model for tumor angiogenesis has recently been proposed and developed in the series of papers [2],[3],[4]. The model consists of a coupled system of ordinary and partial differential equations, which incorporates porous medium diffusion, diffusion by mean curvature and diffusion by reinforced random walk. It consists of two sets of equations, one for the dynamics inside a capillary wall, which is considered as a 1-dimension medium, and one for the dynamics in the adjoining extracellular matrix (abbreviated as ECM), which is viewed as a 2-dimensional medium separating the tumor from the capillary.

In general, the ECM is made of structural proteins. In the model the structural protein is assumed to be of one type, namely fibronectin. A capillary is composed of endothelial cells. In order for the outgrowth of a new capillary to occur, these structural proteins must be broken down. Thus endothelial cells must not only proliferate but also secrete a proteolytic enzyme that degrades the structural proteins surrounding the capillary, thereby providing space for new cells to enter. Continued growth through the ECM requires the same two processes and is directed by chemotaxis, the cellular sensitivity to a chemical or protein gradient.

In modeling this process, concentrations of the tumor emitted growth factor, the structural protein and the degrading proteolytic enzyme are used along with the density of the endothelial cells. The concentrations and densities are functions defined on two different domains, corresponding to the capillary and the ECM. The state of the capillary is given in terms of four functions: v(x,t), the concentration of the angiogenic growth factor; c(x,t), the protease concentration, f(x,t); the fibronectin concentration; and n(x,t), the endothelial cell density. The state of the ECM is also given in terms of four functions, V(x,y,t), C(x,y,t), F(x,y,t) and N(x,y,t), having meanings analogous to their lower case counterparts.

In this paper an algorithm for numerical simulations with the model is described. The simulation results that were reported in \cite{LSNHP} were performed using a MATLAB code in which spatial derivatives were approximated to give a semidiscrete problem that was then integrated using one of the built-in ODE solvers. However, the MATLAB code had several drawbacks, most notably speed of execution. In developing a Fortran 90 code it became clear that one of the reasons for the long excution times is that the problem is highly convection dominated, and cannot be effectively treated as a semidiscrete ODE. The main part of this paper, section 3, describes a method based on characteristics, which is mass conserving, that effectively deals with this feature of the problem. This method is incorpated into an adaptive time stepping algorithm that can only briefly be described in section 4.

References

1. J.Folkman, The vascularization of tumors, Scientic American, 234 (1976), 58-64.
2. H.A. Levine, B.D. Sleeman, M. Nilsen-Hamilton, A mathematical model for the roles of pericytes and macrophages in the onset of angiogenesis: I. The role of protease inhibitors in preventing angiogenesis, Math. Biosci. 168 (2000), 77-115.
3. H.A. Levine, B.D. Sleeman, M. Nilsen-Hamilton, Mathematical modeling of the onset of capillary formation initiating angiogenesis, J. Math. Biol., 42 (2001), 195-238.
4. H.A. Levine, S. Pamuk, B.D. Sleeman, M. Nilsen-Hamilton, Mathematical modeling of capillary formation and development in tumor angiogeneesis: Penetration into the stroma, Bull. Math. Biol., 63 (2001), 801-863.

A Finite Volume Lagrangian Method

A finite volume approach can be used to obtain approximations for most of the differential equations in the modeling system. In our approach the ECM is subdivided into control volumes (i.e. rectangles R_ij), and the equations are integrated over space-time volumes, R_ij x [t_n,t_n+1], to obtain integral identities which are then approximated. Implicit methods are used thoughout since the system is stiff due to widely varying diffusion and chemical rate constants. An alternating direction scheme has been used in approximating the evolution of V(x,y,t) in the ECM. The diffusion of fibronectin is several orders of magnitude small than the other rates, so that the equation for the evolution of the fibronectin, F(x,y,t), can be handle much like an ODE. The paper contains more details on this.

The essential difficulty in obtaining accurate and stable numerical approximations of the total system lies in dealing with the strongly convection dominated diffusion equations for the endothelial cell density in the capillary and in the ECM. Due to the hyperbolic character of the evolution equations for these quantities, a method based on characteristics is used in approximating these equations. This method is the main subject of the paper. Methods of this type have been proposed by several authors for transport problems. However, some approaches fail to produce mass conserving schemes. Our approach, which does yield a mass conserving scheme, is similar to those described in references below. These papers also provide extensive references to the literature.

• T. Arbogast and M. Wheeler, A characteristics-mixed finite element method for advection-dominated transport problems, SIAM J. Numer. Anal. 32 (1995), 404-424.
• J. Douglas, F. Pereira, L. Yeh, A locally conservative Eulerian-Lagrangian numerical method and its application to nonlinear transport in porous media, Computational Geosciences 4 (2000), 1-40.

Examples of simulations performed using this method are included below. The examples consist of a series of graphs giving snapshots in time of the density and concentration functions at various times during a simulation run. The graphs appear over non-dimensionalize domains, the unit interval in the case of the capillary functions and the unit square in the case of the ECM functions. The physical dimensions are L_x = 100 microns and L_y = 20 microns. The number of subdivisions in the x and y directions are M_x = 81 and M_y = 40 respectively. A variable time stepping strategy was used in the simulation run. Initially it was chosen to be dt = 0.01, so that 100 time steps corresponds to 1 hour of model simulation time. As sharp fronts develop the time step must necessarily be changed.

An Example with Protease & Fibronectin Chemotaxis

A Graph at simulation time t = 0.1 hour

The first graph shows the initial influx of growth factor from the tumor side of the ECM (i.e. y = 1). All functions have been started with initial constant states. The ECM functions are initially zero, and the capillary functions are initially: v(x,0) = 0, c(x,0) = 0, f(x,0) = 1, n(x,0) = 1. At this point in time only the growth factor concentration V(x,y,t) has become non-constant. The growth factor from the tumor drives the system as can be seen in the subsequent graphs.

Graphs at simulation time t = 1 hour

Graphs at simulation time t = 2.8 hours

The next set of graphs shows the state of the system after 2.8 hours of model simulation time. The fibronectin surrounding the capillary wall has now been sufficiently degraded for capillary outgrowth to occur. As the new capillary grows accross the ECM it consumes growth factor V(x,y,t) and produces protease C(x,y,t), which in turn degrades the ECM fibronectin F(x,y,t).

capillary functions at t = 2.8 hrs	C(x,y,t) at t = 2.8 hrs
F(x,y,t) at t = 2.8 hrs	N(x,y,t) at t = 2.8 hrs

Graphs at simulation time t = 4.1 hours

The next set of graphs shows the state of the system after 4.1 hours of model simulation time. The outgrowth of the new capillary, represented by the density function N(x,y,t) of the endothelial cells, has continued farther into the ECM region. The cell proliferation and migration is concentrated at the tip of the region growing towards the tumor side of the ECM. The production of protease is seen to be largest at the tip, resulting in the degradation of the ECM fibronectin. This has resulted in a channel being formed. In this simulation run, the chemotactic sensitivity to the ECM fibronectin has be chosen so that a concentration of F(x,y,t) = 0.5 is favorable. Thus endothelial cell movement is directed away from the extremes F(x,y,t) = 0 or 1 and towards F(x,y,t) = 0.5. This is a different choice than the one used in the simulations reported in paper [4] listed above. The graph of the endothelial cell density N(x,y,t) shows the cells tend to line the edge of the channel and form the walls of a new capillary.

C(x,y,t) at t = 4.1 hrs	V(x,y,t) at t = 4.1 hrs
F(x,y,t) at t = 4.1 hrs	N(x,y,t) at t = 4.1 hrs

Graphs at simulation time t = 4.9 hours

These graphs show the endothelial cell density N(x,y,t) after approximately 4.9 hours. The propagating front has a jagged appearance in this graph but has actually been well resolved by the numerics. The jagged appearance is an artifact of the MATLAB graphics package that has been used to produce all of the graphs presented here. The difficulty in plotting is due to the sharp edges of the density surface. The sharpness of the edges of the density surface can be seen by taking traces in vertical planes, x = x_i and y = y_j.

N(x,y,t) at t = 4.9 hrs	N(x,yj,t), j=1,2, at t = 4.9 hrs
	N(xi,y,t), i=1,2, at t = 4.9 hrs

Graphs at simulation time t = 5.6 hours

The last graphs show the endothelial cell density as it approaches the tumor side of the ECM. As seen at earlier times, the cells line the walls of the outgrowing capillary. Blood would be transported through the channel between the walls. This channel is more clearly seen in the rotated graph. Again the sharp edges to the endothelial cell density distribution are apparent. A different view of this is provided by a contour plot.

N(x,y,t) at t = 5.6 hrs	Same Graph Rotated
	A Contour Plot of N(x,y,t) at t = 5.6 hrs

An Example with Growth Factor & Fibronectin Chemotaxis

A Graph at simulation time t = 0.1 hour

In this example the growth factor is confined to a diffusable channel. The first graph shows the initial influx of growth factor from the tumor side of the ECM (i.e. y = 1). The outline of part of the channel can be seen as the edge of the support of the growth factor density function. The channel has two branches on the tumor side of the ECM which merge into one branch as it tranverses the ECM to reach the capillary side. As in the first example, all other functions have been started with initial constant states. The ECM functions are initially zero, and the capillary functions are initially: v(x,0) = 0, c(x,0) = 0, f(x,0) = 1, n(x,0) = 1. At this point in time only the growth factor concentration V(x,y,t) has become non-constant. Again he growth factor from the tumor drives the system.

Graphs at simulation time t = 1 hour

Graphs at simulation time t = 7.3 hours

The next set of graphs shows the state of the system after approximately 7.3 hours of model simulation time. The overall process has been slowed since the diffusion of the growth factor has been hampered by being confined to the channel. The fibronectin surrounding the capillary wall has now been sufficiently degraded for capillary outgrowth to occur. As the new capillary grows accross the ECM it consumes growth factor V(x,y,t) and produces protease C(x,y,t), which in turn degrades the ECM fibronectin F(x,y,t).

V(x,y,t) at t = 7.3 hrs	C(x,y,t) at t = 7.3 hrs
F(x,y,t) at t = 7.3 hrs	N(x,y,t) at t = 7.3 hrs

Graphs at simulation time t = 8.6 hours

The next set of graphs shows the state of the system after approximately 8.6 hours of model simulation time. The outgrowth of the new capillary, represented by the density function N(x,y,t) of the endothelial cells, has reached the location where the channel forks into two branches. The chemotactic senitivity of the endothelial cells to the growth factor can be clearly seen at this point in a contour plot of level curves of the growth factor concentration V(x,y,t) superimposed with those of the endothelial cell density N(x,y,t). As in the first example, the chemotactic sensitivity to the fibronectin has be chosen so that a concentration of F(x,y,t) = 0.5 is favorable. Thus endothelial cell movement is directed away from the extremes F(x,y,t) = 0 or 1 and towards F(x,y,t) = 0.5.

C(x,y,t) at t = 8.6 hrs	V(x,y,t) at t = 8.6 hrs
F(x,y,t) at t = 8.6 hrs	N(x,y,t) at t = 8.6 hrs

Graphs at simulation time t = 9.8 hours

The next set of graphs shows the state of the system after approximately 9.8 hours of model simulation time. The outgrowth of the new capillary has now branched into two outgrowths. Again the chemotactic senitivity of the endothelial cells to the growth factor can be clearly seen at this point in a contour plot of level curves of the growth factor concentration V(x,y,t) superimposed with those of the endothelial cell density N(x,y,t).

Maintain by M.W. Smiley, last updated: April 18, 2003