Study on the tumor-induced angiogenesis using mathematical models.

We studied angiogenesis using mathematical models describing the dynamics of tip cells. We reviewed the basic ideas of angiogenesis models and its numerical simulation technique to produce realistic computer graphics images of sprouting angiogenesis. We examined the classical model of Anderson-Chaplain using fundamental concepts of mass transport and chemical reaction with ECM degradation included. We then constructed two types of numerical schemes, model-faithful and model-driven ones, where new techniques of numerical simulation are introduced, such as transient probability, particle velocity, and Boolean variables.

INTRODUCTION

Angiogenesis, the formation of new blood vessels from pre‐existing vessels, is a vital component of many growth processes,1 including embryogenesis, retinal vasculature, wound healing, tumor growth and numerous vascular diseases. Hypoxic cells release angiogenic factors such as vascular endothelial growth factor (VEGF) to promote angiogenesis. When VEGF reaches nearby pre‐existing vessels, the endothelial cell is induced and activated, leading to the formation of a tip cell with increased motility. Tip cell migrates in response to the gradient of VEGF, followed by migration of the stalk cell behind the tip cell, resulting in the formation of a new blood vessel sprout to carry oxygen and nutrients to hypoxic cells (Figure 1). Sprouting angiogenesis is not the only method to develop new blood vessels. Vasculogenesis, the self‐development of vascular cells, is the other way to promote new blood vessel formation. In the present article, we review sprouting angiogenesis only. For a biological review of vasculogenesis, see Takakura.2

Figure 1

Illustration of angiogenesis. Vascular endothelial growth factor (VEGF) is released by hypoxia cells, diffusing through the ECM, and finally reaching the nearby pre‐existing blood vessel. Tip cells are formed by endothelial cell activation, which is caused by VEGF induction to the pre‐existing vessel. Then, tip cells migrate in response to VEGF gradient chemotactically. Tip cells degrade ECM (orange lines) during this migration. Branching and anastomosis occurs during the angiogenesis process

Some of the earlier angiogenesis models and simulations described the continuous density of network vessels. This type of model lacks details in the formed network of a sprout such as tip branching and anastomosis. An example can be taken from the work of Aubert et al.3 They developed a partial differential system of retinal vascularization that involves the interaction of tip cell, stalk cell, and VEGF. Numerical simulation was carried out on a one‐dimensional (1D) domain. Figure 2 explains how the density of tip cell, stalk cell, and VEGF are defined in a 1D model and simulation. Tip cell density is defined as the cell density of the sprout tip whereas stalk cell density is defined as the cell density of sprouts other than tip. Tip branching and anastomosis are defined as kinetic functions with the form of proliferation and decay terms. By this definition, we cannot see the difference between cell proliferation and tip branching. Tip cell density increases with tip branching applied. The rate of increase depends on the choice of branching rate. Another important phenomenon, anastomosis, is treated similar to tip branching. By defining anastomosis as tip cell decay and stalk cell proliferation, we can see only the increasing or decreasing tip cell and stalk cell density. We do not see the detailed structure of sprouting angiogenesis from their model and simulation. Another example of a continuum model and simulation was performed in two dimensions (2D) by Anderson and Chaplain.4 They involved tip cell, VEGF, and ECM in the model. Tip cell migrates in response to VEGF and fibronectin gradients. We call this chemotaxis and haptotaxis, respectively. They did not consider tip branching and anastomosis in their continuum model. The simulation results show the importance of chemotaxis and haptotaxis in sprouting angiogenesis. Chemotaxis drives tip cell migration directly towards the tumor while haptotaxis creates dispersion of tip cell density.

Figure 2

Density of tip cell, stalk cell, and vascular endothelial growth factor (VEGF) on 1D model are defined as the average density on the green dash‐line for each position

The discrete model is often combined with the continuous model to obtain the detailed structure of a vessel network. Anderson and Chaplain4 carried out this technique to track individual tip cell migration. The tracked path creates a continuous line which is then defined as sprout formation. The tip cell position is updated probabilistically depending on the set of environmental factors (chemotaxis and haptotaxis) at every time step. Another discrete technique was carried out by Bentley et al.5 They introduced an agent‐based model for simulating the tip cell selection process involving tip cell induction which is mediated by VEGF and delta‐like 4 (DII4)‐Notch signaling. Simulation was carried out on a three dimensions (3D) gridded lattice. In the simulation, they defined endothelial cell as ‘EC agents’ which is composed of a small agent (‘memAgent’) on the cell periphery. This memAgent contains receptors of proteins that can promote filopodia extension. Thus, EC agents are capable of filopodia extension and retraction depending on the level of protein surrounding that agent. This filopodia extension formed a continuous line as a new blood vessel. Plank and Sleeman6 developed a lattice‐based model of tumor‐induced angiogenesis. They defined the probability of tip cell migration in each direction as in Othmer and Stevens.7 A reviewed article about recent computational models of sprouting angiogenesis was written by Heck et al.8 They categorized the models into three groups: tip cell migration‐based models; tip cell and stalk cell‐based models, and cell shape dynamics models. The first two groups focus on cell migration based on chemotaxis and haptotaxis. The last group completes the system by involving cell shape dynamics during sprouting angiogenesis. The Cellular Potts model can be used to capture this phenomenon. This is a lattice‐based model that allows simulation of cell behavior based on total energy. Cell behavior such as cell size and shape are translated into an energy equation to become decision points for updating the cell surface at every time step of numerical simulation. Bauer et al.9 developed the first model and simulation of sprouting angiogenesis with cell shape dynamics. This involved a discrete process for migration of cells and a continuous method for the environmental factor (VEGF).

Mathematical study on sprouting angiogenesis, however, integrates biological insights as mathematical formulas to understand complicated events in a clear and quantitative way, to predict morbid states for proposing better therapeutic strategies, and to develop effective tools to create new drugs. This approach was initiated by Anderson and Chaplain4 for tumor‐induced angiogenesis, and was then extended to wound healing10 as well as to retinal vasculature.3, 11, 12 The study is still under progress collaborating with recent developments in cell biology; for example, interactions of tip, stalk, and mural cells under the control of cellular molecules such as Angiopoietin 1, Angiopoietin 2, platelet‐derived growth factor B (PDGF‐B), and so forth.13

The purpose of the present review is to describe their basic ideas, expecting the reader to see the benefits of these studies, recognizing the core part of mathematics used in the models. The reader is expected to also become a good user or collaborator and, furthermore, one of the experts engaged in this new field of science, mathematical oncology. To this end, we pick up the classical model,4 modify it based on later biological knowledge, and propose new numerical schemes which lead to realistic computer graphics images as in http://bodytokyo.jp/. We review the comparison of numerical simulation results with experimental data to show the relevance of the model.

BIOLOGICAL INSIGHTS

We begin with a review of several biological insights recently established on tumor‐induced angiogenesis. First, three types of VEGF are noticed; that is, VEGF120, VEGF165, and VEGF189. VEGF120 and VEGF189 are binding and unbinding ECM, respectively, and VEGF165 takes an intermediate profile. Strong expressions of VEGF120 and VEGF165 are observed in angiogenic tissues.14, 15 ECM‐binding VEGF is released with ECM degradation by MMP, which is thought to be the origin of the VEGF gradient.16, 17, 18 In total, chemotaxis acting on the tip cell is caused by VEGF gradient, and this gradient is due to ECM degradation by MMP.

Second, remodeling of ECM by tip cells19, 20, 21 and also haptotaxis induced by ECM degradation by MMP22 are observed. Finally, there is MMP upregulation inside the endothelial cell activated by VEGF fragments,18, 23, 24, 25 particularly inside the tip cell.22

FUNDAMENTAL CONCEPTS OF MATHEMATICAL MODELING

Here we describe two basic concepts in mathematical modeling. The first is ordinary differential equation (ODE) describing the amount balance. Supply and consumption with the rate α of the quantity u = u(t) varying with the time t, are described byrespectively, where , while

respectively, indicate growth and decay u with rate β. The quadratic non‐linearity, furthermore, is used to describe the interaction of two types of particles, as in the chemical reaction inside the solution, A + B ⇀ P(k); that is,

This model leads to mass conservation, and, with the constants a and b determined by initial values, which is used to simplify the reaction network in a quasi‐stationary state. For example, if the receptor‐ligand processis in equilibrium on the plasma membrane, the formula of Michaelis‐Menten holds aswhere c  =  [RW] + [R] is a constant determined by the initial value and γ = l/k is the equilibrium constant.

The second concept is the gradient. Given the scalar field φ, its gradient ▿ represents the vector field with the direction maximizing the growth of φ with the length of its rate. The gradient operator ▿ is thus defined, which leads to the divergence ▿j of the vector field j. If ρ = ρ(x, t) represents the mass density varying in t, its flux j = j(x,t) is defined bywhere ω is an arbitrary domain, v is the outer normal unit vector, and dS is the surface element. Then the divergence formula of Gauss induces the equation of conservation,

In the diffusion equationtherefore, the flux j of ρ is proportional to ‐ ▿ρ. Here and henceforth, d > 0 represents the physical constant associated with the material. Some other positive constants are denoted by α, β, γ, δ, and so forth, as before. In the transport theory, the variation of ρ in time is described in accordance with the velocity of particles, denoted by v = v (x, t) and, generally, the equation in the form ofis called the transport equation, wheredenotes the material derivative. Under this flow, the position x = x (t) of the particle is subject toand if ω denotes the region made by the particles at time t which are in ω at t = 0, Liouvilles's theorem implies

Mass conservation is thus reduced toand the relation j = ρv follows, which indicates that the flux of particle density is the product of itself and its velocity.

TIP CELL MODEL

Here we modify the classical model4 using recent insights of cell biology and mathematical modeling. The modification is based on recent biological knowledge reviewed in the section Biological Insights. First, tip cell is assumed to be continuously distributed and hence is represented by n = n (x,t) with x and t representing space and time variables, respectively. It is subject to diffusion and also chemotaxis by VEGF and haptotaxis caused by ECM degradation. This process is represented bywhere c = c (x, t) = c(x, t) and are VEGF concentration and ECM (fibronectin) density, respectively, and v represents the velocity of n‐particles other than the diffusion. The second term on the right‐hand side of v describes one of the driving forces of tip cell movement, called haptotaxis, the invasion to ECM of the tip cell. The first term on the right‐hand side of v , however, represents chemotaxis, so that indicates chemotactic sensitivity. For the single cell the formdue to is standard,26, 27 (page 262) but may be represented as a constant at the tissue level.

ECM density change is under the control of the tip cell, which is two‐fold; that is, ECM degradation and remodeling. This process is formulated byusing the material derivative with constant rate μ on the left‐hand side of the first equation, where v denotes the velocity of n defined by the second equation. The first term on the right‐hand side of the first equation represents the remodeling of ECM by the tip cell, whereas the second term is concerned with ECM degradation by the tip cell. This process is activated by the signaling caused by the VEGF fragment c * = c 0 − c in accordance with MMP proliferation inside the cell, where is the initial distribution of ECM. Hence, in the above model of ECM degradation, MMP concentration is replaced by tip density. Here, it is reasonable to assume η to bedue to the ligand‐receptor binding process, using the equilibrium constant γ, which, however, may be a constant at the tissue level.

The localized VEGF attached to ECM is thought to be the origin of its gradient. This type of VEGF is destroyed by MMP which is identified with the tip cell as above. Hence, it follows that

Totally, only n is involved by the diffusion in this system, which should be provided with the boundary condition. It is reasonable to assume the null‐flux condition here, but may be replaced by the Neumann zero condition as c and f can remain as constants near the boundary in a short time because they are subject to hyperbolic equations. Mathematical justification of this system is similar to several models associated with diffusion and chemotaxis.28, 29

MODEL FAITHFUL DISCRETIZATION

Here we describe the method of finite difference for the simple case of 1D‐interval, I = [0,1], using the uniform mesh h defined by hN = 1, where N is a large integer. Usually, is taken as the main lattice, where even and odd extensions are taken at the end points x = 0, 1 for Neumann and Dirichlet boundary conditions, respectively. On this lattice we have two ways of approximation of the first derivative,

which yields the center difference approximation of

One may take the mean

for approximation of the first derivative on the main lattice. The other way valid to the convection equation

is the upstream difference,with , where i and k represent the lattice and the time step indices, respectively. The mean for approximation of the first derivative may be replaced by

using the sub‐lattice denoted by

To solve the equation of conservation

numerically, we use the non‐uniform time mesh combined with the uniform space mesh,where denotes the space discretization obtained by the above methods and 0 < θ < 1 is a constant, so that the mixed Euler difference scheme is applied for time discretization. Writing the above scheme simply likewe apply Varga's theorem to A .30 It says that any irreducible, diagonally dominant square matrix A, of which diagonal entries are positive and the others are non‐positive, is reversible with A −1 composed of only positive entries. The assumption of this theorem holds to A = A if τ is sufficiently small, with the bound calculated by the quantities at the k‐step. By this adaptive mesh, the positivity and the total mass of n are maintained step by step, provided by the null‐flux boundary condition.31, 32 We obtain a similar fact foradapting the center difference scheme for the diffusion term. Once this mechanism of numerical stabilization is confirmed, even forward Euler scheme θ = 0 is effcient for actual numerical simulations.33

SCHEME FOR HYBRID SIMULATION

Writing the above scheme as

we can regard as the transient probabilities of the particle on the sites to that i at each step, because it follows that for and from the positivity preserving and the total mass conservation. This scheme is easily extended to the case of two‐dimensional space.

Introducing the above transient probabilities is a fundamental concept for hybrid simulation described below. Here, we formulate discrete total velocity. First, the velocity v of n other than the diffusion is determined by c and, which is regarded as the environment variables. These variables are subject to the first‐order equation with the material derivative

Then, the first term on the above right‐hand side is discretized by the forward Euler scheme to proceed to the next time scheme, while upwind finite difference is used for the second term. For this purpose, we define the velocity v by the formula

where . In the case of two‐dimensional space, the right‐hand side is composed of four terms involving the unit vectors . In fact, we recall that h and τ represent space mesh and time mesh, respectively, and the other term on the right‐hand side indicates the dimensionless quantity indicating the position of n at the next time step in probability.

The other concept of our discretization for hybrid simulation is the use of the Boolean variable to tip cell density, introduced by.11 Thus, the variable n is localized at several lattice points, denoted by n*, and we replace n by n* in the right‐hand side on the transport equations concerning the environment variables; that is, c and. We thus end up with a closed full discrete system concerning n*, c and.

NUMERICAL SIMULATION

In this section, we describe the numerical technique of hybrid simulation and show examples of simulation results. We adopted the parameter values from Anderson and Chaplain.4 The domain is square with length 1 (non‐dimensional). We use d  = 3.5 × 10−4, χ = 0.38, α = 0.05, γ = 0.1, and δ = 0.1.4 We assume β = 1, μ = 10, and μ = 10 to satisfy the stability condition of the numerical scheme. For this paper, we take χ to be constant. Taking χ to be dependent of VEGF concentration leads to decreasing of chemotactic sensitivity. Tip cell will become less sensitive to high VEGF, leading to deceleration of tip cell migration. Detail of simulation using VEGF dependent‐χ can be seen in Minerva.33

In the simulation setting, tip cell is assigned as a single point. We set five tip cells in the beginning which are denoted as five red dots at the left boundary (Figure 3A). We choose initial distribution of ECM to be uniformly distributed,(Figure 3B) and VEGF to be gradually distributed with the form(Figure 3C).

Figure 3

Initial condition setting of (A) five tip cells position, B, ECM distribution, and C, vascular endothelial growth factor (VEGF) distribution

At every time step, we decide the next movement of each tip cell by using transient probability written similar to Eqn (29). The choices of tip cell movement are to stay, to move to the left, to the right, down, or up. Formation of blood vessels is visualized by tracing tip cell movement in computer graphics. We also put in the branching and anastomosis rules. In the former, each tip cell is subject to branching with the probability subject to VEGF density denoted by c.4 We add a new tip cell with a random choice of direction if branching occurs (Figure 4B). During the simulation, a tip cell can meet another tip cell or formed vessel to form loop formation. This is what we define as anastomosis in numerical process. In the case of a tip cell that meets another tip cell, we select one tip cell randomly to maintain the model‐driven scheme at the next time step. We remove a tip cell if that tip cell meets the formed vessel (see Figure 4A for detail). This definition of branching and anastomosis can also be applied to the 3D domain of simulation. More choices of direction for tip cell migration give room for the tip cell to branch. Remember that the rule of branching influences the occurrence of tip branching. The choice of initial distribution of ECM and VEGF, and position of tip cells at initial of simulation also affect this two phenomenon. An example of 3D simulation was carried out by McDougall et al.11 They successfully discovered branching and anastomosis during the 3D simulation.

Figure 4

A, Types of anastomosis. (Left) Tip cell meets other vessel or what we call tip to sprout anastomosis. (Right) One tip cell meets another tip cell or what we call tip to tip anastomosis. B, Branching form possibilities

At the same time step as in tip cell migration, distribution of ECM and VEGF is also updated by solving Eqn (16) and Eqn (18) numerically using Euler approximation. Velocity v is calculated using Eqn (31). Tip cell n is valued as 1 or 0 in Eqn (16) and Eqn (18), depending on the existence of tip cell (1 is for presence).

We now show some examples of simulation results. Figures 5 and 6 are snapshots of numerical simulation. Red color dot represents tip cell existence. The black line is the tracked path created by tip cell migration. We define this as a newly formed vessel. By this time, a vessel is formed and approaching the right boundary where high VEGF occurs (Figure 5A). Branching and anastomosis occurs after the vessels reach position x = 0.2. This two phenomenon can occur depending on the choice of VEGF density and branching probability. Placing a different profile of VEGF and ECM affects the occurrence of branching and anastomosis.33 Furthermore, ECM degradation and remodeling are observed during tip cell migration in Figure 5B and 6. The color changes on the area of a newly formed vessel are detected as this two phenomenon. This creates an ECM gradient difference for the tip cell to migrate. Minerva33 carried out alterations to each parameter in the numerical simulation. She also used two different density profiles of VEGF and ECM to further examine the difference of vessel growth on the alterations. In the absence of ECM remodeling and ECM degradation, vessel growth has no significant changes. Using a different profile of ECM distribution at the initial step, vessel growth can be inhibited. VEGF induction on the pre‐existing vessel results in fibronectin leakage from the vessel. This fibronectin bound to ECM forms an ECM gradient where the high concentration is on the pre‐existing vessel area. With chemotaxis term applied in this ECM gradient, tip cell migrates forward and backward, creating a ‘brush border’, and never reaching the right boundary where the tumor is placed in the model of tumor‐induced angiogenesis.4 Sholley et al.34 observed that vascularization in the tumor did not occur in the absence of cell proliferation. Based on this fact, Anderson and Chaplain added cell proliferation to the simulation process. As a result, the vessel now reached the right boundary, penetrated the tumor where we assume it to be placed.4

Figure 5

Snapshots of vessel growth on (A) vascular endothelial growth factor (VEGF) distribution and B, ECM distribution in time

Figure 6

Snapshots of ECM distribution in time

Alteration of VEGF distribution and equation is more interesting. Placing the source of VEGF, which is regarded as tumor, at the center of the right boundary gives a different VEGF gradient profile from Figure 2 (right). The pattern of vessel networks follows the pattern of VEGF gradient.33 In the absence of VEGF degradation, setting δ = 0, vessel growth is decreasing. The vessels stop growing once they reach the highest density of VEGF.33 This may result in the failure of complete vascularization inside the tumor.

COMPARISON WITH EXPERIMENTAL DATA

There are many earlier studies on mathematical models of sprouting angiogenesis involving experimental data as comparison. Aubert et al.3 used retinal vasculature data in mice as comparison for simulation of the capillary tip and astrocytes migration on 1D. They detected the position of the capillary and astrocyte tip at every time increment and compared them with experimental data. Their data match very well. McDougall et al.11 further compared their 2D model and simulation results of retinal vasculature with experimental data in mice. Not only tip position data, they successfully compared the structure of a vessel network with a real figure. The radii changes that are related to their angiogenesis system creates heterogeneity of vessel radii which makes visualization of computer graphics more realistic.

CONCLUSION

Mathematical modeling of sprouting angiogenesis at the tissue level using partial differential equations and visualization trials using hybrid simulation with transient probabilities and Boolean variables have been computed. However, these are based on empirical rules and measured values. In this article, we established a strategy to reflect the remodeling and recent biological knowledge that are directly consistent with numerical simulation results. We proposed a modification of the classical model4 on tumor‐induced angiogenesis based on recent developments in biological and mathematical theories. We formulated ECM degradation and penetration by tip cell using the transport theory. We clarified theoretically that the elements missed in the simulations, which included alterations in the model, that have been done so far do not become practical obstacles because no specific difference was found in numerical simulation results. Model‐faithful and model‐driven discretized schemes were introduced, where particle velocity is formulated besides the Boolean environment variables in accordance with the transient probabilities of the particle variable. Using our discretization scheme that guarantees positivity preservation and mass conservation, we showed that the transient probability adopted in hybrid simulation is necessarily valued between 0 and 1. Deterministic and stochastic hybrid simulations then visualized realistic movements of tip cells.

CONFLICT OF INTEREST

Authors declare no conflicts of interest for this article.