Ensemble Averages , Soliton Dynamics and Influence of Haptotaxis in a Model of Tumor-Induced Angiogenesis

In this work, we present a numerical study of the influence of matrix degrading enzyme (MDE) dynamics and haptotaxis on the development of vessel networks in tumor-induced angiogenesis. Avascular tumors produce growth factors that induce nearby blood vessels to emit sprouts formed by endothelial cells. These capillary sprouts advance toward the tumor by chemotaxis (gradients of growth factor) and haptotaxis (adhesion to the tissue matrix outside blood vessels). The motion of the capillaries in this constrained space is modelled by stochastic processes (Langevin equations, branching and merging of sprouts) coupled to continuum equations for concentrations of involved substances. There is a complementary deterministic description in terms of the density of actively moving tips of vessel sprouts. The latter forms a stable soliton-like wave whose motion is influenced by the different taxis mechanisms. We show the delaying effect of haptotaxis on the advance of the angiogenic vessel network by direct numerical simulations of the stochastic process and by a study of the soliton motion.


Introduction
When an avascular tumor surpasses some 2 mm, normal tissue vasculature can no longer support its growth.The tumor cells then lack oxygen and nutrients and produce vessel endothelial growth factors and other tumor angiogenic factors (TAF).These substances reach nearby blood vessels that emit blood vessel sprouts as a response.The vessels reach the tumor bringing oxygen and nutrients that allow it to get bigger and proliferate.In this case, a tumor induces the growth of blood vessels, a process called angiogenesis [1].While angiogenesis is a natural process responsible for organ growth and repair, imbalance between stimulating and inhibiting factors contributes to numerous malignant, inflammatory, ischaemic, infectious, and immune disorders [2].While tumor cells do not become malignant due to angiogenesis, this process promotes tumor progression and metastasis [1,2].Angiogenesis encompasses multiple scales, from cellular and subcellular micron scales to the mesoscopic dynamics of millimeter sized blood vessels.Involved time scales also vary widely from seconds to days [3,4].Many aspects of angiogenesis have been discovered by combining laboratory experiments and results from computational models; see the review paper [4].
The complexity of angiogenesis requires posing hybrid models that combine a detailed description of individual vessel sprouts (or even individual cells) with continuum descriptions of the growth factors and of the extracellular matrix (ECM) over which vessels advance in a constrained space.Most work consists of comparing numerical simulations of these hybrid models to experimental results, whereas deeper analyses are lacking.Among hybrid models, the simplest ones are called tip cell models [4].They assume that the endothelial cells (ECs) at the tip of a blood vessel are highly motile, do not proliferate and respond to the growth factors and chemical cues that lead them to the tumor.Other endothelial (stalk) cells proliferate and doggedly follow tip cells, constructing the vessel in the meantime on their wake.We can treat the tips of the sprouting blood vessels as particles subject to different forces (chemotaxis, haptotaxis, . . . ) that link them to continuum fields (representing growth factors, inhibitors, ECM, . . . ) and consider their trajectories as the expanding angiogenic vessel network.It turns out that vessel tips may branch (thereby creating new active vessel tips) or merge with existing vessels when they reach them, in a process called anastomosis.In the latter case, the tip cells of the merged vessel cease to be active.
The earlier stochastic tip cell model was proposed by Stokes and Lauffenburger in 1991.They considered the capillary sprouts as particles of unit mass subject to chemotactic, friction and white noise forces [5,6].Associated to each sprout, its cell density satisfies a rate equation that takes into account proliferation, elongation, redistribution of cells from the parent vessel, branching and anastomosis.They did not consider the depletion effect that advancing sprouts would have on the TAF concentration.Later tip cell models combined a continuum description of fields influencing cell motion (chemotaxis, haptotaxis, . . . ) with random walk motion of individual sprouts that experience branching and anastomosis [4].Capasso and Morale [7] used ideas from these approaches to propose a hybrid model of Langevin-Ito stochastic equations for the sprouts undergoing chemotaxis, haptotaxis, branching and anastomosis coupled to reaction-diffusion equations for the continuum fields.In this model, the evolution of the continuum fields is influenced by the growing capillary network through smoothed (or mollified) versions thereof [8].Capasso and Morale also attempted to derive a continuum equation for the density of moving tip cells from the stochastic equations but could not account for branching and anastomosis [7].Recently, we have derived a deterministic description for the average vessel tip density of a related hybrid stochastic tip cell model for tumor-induced angiogenesis that includes chemotaxis but not haptotaxis [9,10].We have shown that the density of active vessel tips advances toward the tumor as a soliton-like wave [11,12].Besides branching and anastomosis, in this model, tip cells experience a chemotactic force that leads them toward increasing values of the TAF concentration, plus friction and Brownian motion.The TAF diffuses and is consumed by the motion of the tip cells.This picture ignores that the ECs have to move over the ECM.They do so by using adhesion (haptotaxis) in response, e.g., to fibronectin gradients and releasing matrix degrading enzymes to corrode the ECM.One simple way to model haptotaxis is to model fibronectin (FN) concentration and matrix degrading enzymes (MDE) as continuum fields satisfying reaction-diffusion equations, as postulated in [7].Space inhomogeneity of the ECM-fibronectin network affects the speed and density of sprout formation, particularly when we consider cellular and sub-cellular length scales.In our model, we only consider larger mesoscopic scales.Sprout branching depends only on the local TAF concentration.As TAF is consumed by the moving sprouts, inhomogeneity induced by the encroaching vessel network is considered albeit in a primitive fashion.Within the bounds of our mesoscopic model, we could include additional heterogeneity effects by making the branching rate dependent on fibronectin and by using an initial fibronectin distribution that mimics a pre-established ECM distribution.We will ignore these considerations for the sake of simplicity and because their influence on soliton motion is small.To properly account for these effects, we would need detailed models that include them at the cellular level and to derive mesoscopic models from them.Such analyses are work for the future.
Other models assume that vessel tips follow a reinforced random walk [13][14][15][16], which, within appropriate limits, yields Langevin equations [17].In more detailed models, endothelial cells change form and move over the extracellular matrix by Monte Carlo dynamics of a cellular Potts model [18].In appropriate long time limits, the corresponding discrete time random walk reduces to Langevin equations [17,19].Thus, the ideas presented in this paper might be useful when considering other related angiogenesis models.
There are other processes that are relevant in later stages of angiogenesis and that we do not consider in the present work.Once a vascular network is being created, blood flows through the capillaries, anastomosis enhances the flow in some of them and secondary angiogenesis may start in new vessels.Pries and coworkers have modeled blood flow in a vascular network and the response thereof to changing conditions such as pressure differences and wall stresses [20,21].This response may remodel the vascular network by changing the radii of certain capillaries and altering the distribution of blood flow [20,21].McDougall, Anderson and Chaplain [22] have used this formulation to add secondary branching from new capillaries induced by wall shear stress to the original random walk tip cell model [13].Blood flows according to Poiseuille's law, mass is conserved, there are empirical expressions for blood viscosity and for the wall shear stresses, and radii of capillaries adapt to local conditions.Secondary vessel branching may occur after the new vessel has reached a certain level of maturation and before a basal lamina has formed about it [15,22].During such a time interval, the probability of secondary branching increases with both the local TAF concentration and the magnitude of the shear stress affecting the vessel wall.McDougall et al.'s model can be used to figure out how drugs could be transported through the blood vessels and eventually reach a tumor [15,22].In dense vessel networks, secondary branching may have little effect on the number of active tips at a given time, as anastomosis could eliminate secondary branches quickly.Thus, we may ignore secondary branching when considering the density of active tips in such networks.Of course, we cannot ignore it when describing blood flow and network remodeling.One missing feature of these angiogenesis models that take blood flow into account seems to be pruning.It is known that capillaries with insufficient blood circulation may atrophy and disappear.Pruning such blood vessels is an important mechanism to achieve a hierarchical vascular network with optimized transport, as recent global optimization and adaptation algorithms have shown for a number of biological vascular systems [23,24].
The purpose of this paper is to derive a deterministic description of tumor-induced angiogenesis using the Capasso-Morale stochastic tip cell model that includes both chemotaxis and haptotaxis [7].We do not consider blood circulation and later stages of remodeling and pruning the angiogenic network.We find that the density of active vessel tips obeys an integropartial differential equation coupled to reaction-diffusion equations for the concentrations of TAF, FN and MDE.After an initial stage, the density forms a stable two-dimensional lump (angiton) that moves toward the tumor.A one-dimensional profile of the angiton is a soliton-like wave whose size and velocity change slowly according to differential equations that depend on spatial averages involving the TAF, FN and MDE concentrations and describe the advancing angiogenic network.Thus, haptotaxis can be incorporated to a similar framework as the model including only chemotaxis.Its effects can be ascertained by studying the motion of the soliton as given by the solution of the collective coordinate equations.We confirm this picture by comparing the results to simulations of the stochastic process.

The Stochastic Model
The stochastic part of the hybrid model consists of Langevin equations for the extension of vessels, a tip branching process, and anastomosis (namely, the destruction of tips that merge with existing vessels).For details on the latter two ingredients, we refer to [10].Let X i (t) and v i (t) be the position and the velocity of the ith tip at time t.The extension of this vessel is given by for t > T i (T i is the random birth time of the ith tip).At time T i , the velocity of the newly created tip is selected out of a normal distribution with mean v 0 (|v 0 | is measured in units of 40 µm per hour [5] according to Table 1) and variance σ 2 v , while the probability that a tip branches from one of the existing ones during an infinitesimal time interval (t, t + dt] is proportional to ∑ N(t) i=1 α(C(t, X i (t)))dt, where N(t) is the number of tips at time t, α(C) = α 1 C/(C + C R ), and α 1 and C R are positive constants.The tip i disappears at a later random time Θ i , either by reaching the tumor or by anastomosis, i.e., by meeting another capillary.At time t, anastomosis for the ith tip occurs at a point x such that X i (t) = x and X j (s) = x for another tip that was at x previously, at time s < t.In Equation ( 1), W i (t) are i.i.d.Brownian motions, and k (friction coefficient) and σ are positive parameters (see Table 2).
Table 1.Units to non-dimensionalize the model equations.We assume that tip cell migration is controlled by two main mechanisms, namely, (i) chemotaxis, in response to a TAF, which is released by tumor cells and has concentration C(t, x); and (ii) haptotaxis, in response to the gradient of FN, which is present in the ECM and has concentration f (t, x).The force F appearing in Equation ( 1) is then where D 1 , γ 1 , and D 2 are positive parameters (see below).A similar description can be found in [7].The fields appearing in Equation ( 2) satisfy continuum equations as follows: • Tumor angiogenic factor.The TAF diffuses and is consumed by advancing vessel tips according to [10] Here, d c and η c are positive parameters (see below), while δ σ x (x) is a regularized delta function (e.g., a Gaussian with standard deviation σ x ).The TAF is consumed in the process of enlarging the capillary: at the time interval (t, t + dt), the ith tip advances to v i (t) dt, which adds a length |v i (t)| dt to the capillary.Thus, the consumption term should be proportional to For the slab configuration (see below) and the parameter values considered in this paper, the difference with the smaller sink term included in Equation ( 3) is very small.The sink term is local in space because the region about the tips affecting TAF consumption is of the same order as the tip cell size (microns), which is very small compared to the millimeter lengths described by our model.

•
Fibronectin.Following [7], FN is an adhesive substance attached to the ECM that does not diffuse.In addition, FN is also produced by vessel tips and its degradation depends on the MDE concentration m(t, x) produced by the ECs themselves.Hence, we have where β f and η f are positive parameters (see Section 3).
• Matrix degrading enzyme.As mentioned in [15], TAF and FN bind to specific membrane receptors on ECs.The latter produce a MDE, which enhances their attachment to FN contained in the extracellular matrix.The ECs are then able to exert the traction forces required to propel themselves during migration.Here, we consider the model proposed in [7], which includes diffusion, production, and degradation of MDE as where d m , β m , and η m are positive parameters (see Section 3).A more elaborated model can be found in [25].
Initial and boundary conditions for the TAF field C have been proposed in [10].The initial and boundary conditions for Equations ( 4) and ( 5) will be specified later.
In order to obtain in the next section dimensionless equations to perform our numerical simulations, we consider the values shown in Table 1.The parameters appearing in the equations of the model are taken from [10] and listed in Table 2.
Remark 1.The reference values for x, v, t, and C are those considered in [10], as we will use the same numerical setting.These values are consistent with Stokes and Lauffenburger's experiments [5].

Remark 2.
As suggested at page 866 of [13], we take a reference concentration f R for FN equal to the reference concentration C R for the TAF field.

Deterministic Description for the Average Density of Active Vessel Tips
There is a counterpart to the stochastic model for the densities of vessel tips and the vessel tip flux, defined as ensemble averages over a sufficient number N of replicas (realizations) ω of the stochastic process: As N → ∞, these ensemble averages tend to the tip density p(t, x, v), the marginal tip density p(t, x), and the tip flux j(t, x), respectively.Numerical simulations are based on the Euler-Maruyama method to solve the stochastic Equation ( 6), explicit finite-difference schemes to solve Equations ( 7)-( 9) with appropriate boundary conditions as specified below, and stochastic branching and anastomosis of vessel tips as described above.Initially, there are 2N 0 vessel tips whose location on the primary vessel depends on the geometry of the angiogenic system with appropriately distributed random velocities (see below).We carry out numerical simulations for N different randomly selected initial conditions.The ensemble averages are calculated at each time step by using Equations ( 10)- (12) with N = 400 [10].
In [10], it is shown that the angiogenesis model has a deterministic description for the average density of active vessel tips, p(t, x, v), (in brief, deterministic description) based on the following equation: The two first terms on the right-hand side of Equation ( 13) correspond to vessel tip branching and anastomosis, respectively.
Here, we consider a slab geometry with the primary vessel located at x = 0 and the tumor at x = 1, as illustrated in Figure 1.For the sake of simplicity, we assume that v 0 in this configuration is independent from the direction of the parent vessel and forms a small angle π/10 with the x-axis.Other possible choices change the birth term a little in Equation ( 13), but their effects on the marginal density p(t, x, y) are small for the slab configuration.The boundary conditions for the tip density are [9] p where p + = p for v > 0 and p − = p for v < 0, v = (v, w).An absorbing boundary condition p = 0 on the tumor surface would be more realistic than Equation ( 17).This would be computationally more costly as we would need to include a slab that extends beyond x = 1.However, the difference with the results of Equation ( 17) would be appreciable at the last stage when the vessel tips arrive at the tumor, which we do not study specifically in the present paper.The tip flux density at x = 0, j 0 (t, y), obtained by integrating Equation ( 16) for the vector velocity v 0 = (v 0 , w 0 ), with |v 0 | = 1.Different from [9], we have included the step function θ(τ − t) in Equation (19).With τ = ∞ as in [9,10], the primary vessel keeps injecting tip density for all time.However, this may be artificial, as the primary vessel does not inject any more vessels after t = 0 + in many experiments on early stage angiogenesis.Then, τ in Equation ( 19) may be a small time of the order of the time step used in a numerical code.The original boundary condition in [9] did not include the unit step function θ(τ − t) and, as a consequence, the deterministic description given by the tip density equation and its boundary conditions had an artificial injection of tip density at x = 0 for all t > 0.
x/L y/L The anastomosis coefficient, Γ, has to be fitted using the stochastic description, [10].For τ = ∞, there is an extra tip density injection at the primary vessel as compared to the case of a finite and very small τ > 0. This implies that the anastomosis coefficient in the deterministic description is larger for τ = ∞ (as in [10]) than for the small τ > 0 that resembles more accurately the stochastic description.As the anastomosis coefficient affects the shape of the soliton appearing in the next section, we have selected a smaller value for Γ that is 75% of the deterministic anastomosis coefficient in [10], i.e., Γ = 0.109.Equations ( 7)-( 9) become where j(t, x) is the current density (flux) vector at any point x and any time t ≥ 0, The boundary conditions for the TAF concentration are (a = 1.1 and b = 0.3 is half the tumor width) and lim y→±∞ C = 0.The boundary conditions for MDE are zero-flux at x = 0 and x = 1 and lim y→±∞ m = 0.As initial conditions, we use m(0, x, y) = 1, for c = 1.5, d 2 = 0.45 [10,13].The initial tip density is with l x = 0.06, l y = 0.08 in our nondimensional units.This initial condition corresponds to the following initial condition for the stochastic process: There are 2N 0 = 20 initial tips at x = 0, with vertical positions ±y i equally spaced on the interval [−L y , L y ] (L y = 0.3), whose initial velocities are normally distributed about v 0 with standard deviation 1.Note that, as l x and l y tend to zero, Equation (26) becomes The system of Equations ( 13) and ( 20) with appropriate initial and boundary conditions and without haptotaxis has a unique solution, as proved in [26,27].

Soliton and Collective Coordinates
In the overdamped limit, it is possible to derive the following reduced equation for the marginal tip density [12], For constant F = (F x , F y ), µ, and zero diffusion, 1/β = 0, Equation ( 28) has the following soliton-like solution on −∞ < x < ∞: where K is a constant.Numerical simulations on a slab geometry show that the marginal tip density evolves toward (30) after an initial stage [11,12].The profile (30) ceases to be a good approximation to the solution of (28) when X(t) is close to x = 1.For the rest of this section, we shall consider the intermediate stage after the soliton has formed and before it reaches the tumor at x = 1.
A small diffusion and slowly varying continuum fields C, f , m produce evolution equations for the collective coordinates K, c, and X.Then, the marginal density profile at y = 0 can be reconstructed from (30) with spatially averaged F x and µ [12].Note that ps is a function of ξ = x − X and also of x and t through C(t, x) and f (t, x), We assume that the time and space variations of C and f , which appear when ps is differentiated with respect to t or x, produce terms that are small compared to ∂ ps /∂ξ.As explained in [12], we shall consider that µ(C) is approximately constant, ignore ∂C/∂t because the TAF concentration varies slowly (the dimensionless coefficients κ c and χ c appearing in the TAF Equation ( 20) are very small according to Table 3) and ignore ∂ 2 ps /∂i∂j, where i, j = K, F x .We now insert (30) into ( 28), thereby obtaining (32) Equation ( 28) with 1/β = 0 and constant F and µ, has the soliton solution (30).Using this fact, (32) becomes We now find collective coordinate equations (CCEs) for K and c.As the lump-like angiton moves on the x-axis, we set y = 0 to capture the location of its maximum.On the x-axis, the profile of the angiton is the soliton (30).We first multiply (33) by ∂ ps /∂K and integrate over x.We consider a fully formed soliton far from primary vessel and tumor.As it decays exponentially for |ξ| 1, the soliton is considered to be localized on some finite interval (−L/2, L/2).The coefficients in the soliton formulas (30) and the coefficients in (33) depend on the TAF, FN and MDE concentrations at y = 0; therefore, they are functions of x and time and get integrated over x.TAF, FN and MDE vary slowly on the support of the soliton, and therefore we can approximate the integrals over x by The interval I over which we integrate should be large enough to contain most of the soliton, of extension L. Thus, the CCEs hold only after the initial soliton formation stage.Near the primary vessel and near the tumor, the boundary conditions affect the soliton and we should exclude intervals near them from I. We shall specify the integration interval I in the next section.Acting similarly, we multiply (33) by ∂ ps /∂c and integrate over x.From the two resulting formulas, we then find K and ċ as fractions.The factors 1/L cancel out from their numerators and denominators.As the soliton tails decay exponentially to zero, we can set L → ∞ and obtain the following CCEs [12]: In these equations, all terms varying slowly in space have been averaged over the interval I.The last two terms in (34) are odd in ξ and do not contribute to the integrals in (36) and (37), whereas all other terms in (34) are even in ξ and do contribute.Most integrals appearing in (36) and (37) are calculated in [12].The resulting CCEs are in which the functions of C(t, x, y), m(t, x, y) and f (t, x, y) have been averaged over the interval I after setting y = 0. We expect the CCEs (38) and (39) to describe the mean behavior of the soliton whenever it is far from primary vessel and tumor.

Numerical Results
In this section, we shall solve numerically the full stochastic model and obtain the marginal vessel tip density, the TAF, FN and MDE densities by ensemble averages as explained in [10].Numerical simulations indicate that, for the selected initial configuration of fibronectin that decreases away from the primary vessel, haptotaxis delays the advance of the vessel network compared to the model with only chemotaxis (see Figure 1).This is further illustrated by Figure 2 that compares the profiles of the marginal tip density (at y = 0) including and excluding haptotaxis for two different times.Figure 3 shows how the advancing angiogenic network consumes fibronectin at a smaller rate as it approaches the tumor because the MDE concentration decreases rapidly with time according to (9) and Table 3.
Numerical simulations show that the density of active tips approaches the soliton described in the previous section after an initial formation stage and before the soliton is too close to the tumor at x = 1.The stable soliton instantaneously adapts its shape and velocity according to the solution of the CCEs (38) and (39).From these simulations, we obtain the evolution of the soliton collective coordinates and reconstruct the marginal tip density at y = 0 from (30).As shown in Figure 4, the soliton produces a simple description of tumor induced angiogenesis that agrees with numerical simulations of the stochastic process, provided the maximum of the marginal tip density is not close to the tumor.The simulations show that the soliton is formed after some time t 0 = 0.2 (10 h) following angiogenesis initiation.To find the soliton evolution afterwards, we need to solve the CCEs (38) and (39) whose coefficients are spatial averages (40) that depend on the TAF, FN and MDE concentrations and their derivatives calculated at y = 0.The averaging interval I should exclude regions affected by boundaries.We calculate the spatially averaged coefficients in (38) and (39) by: (i) approximating all differentials by second order finite differences; (ii) setting y = 0; and (iii) averaging the coefficients from x = 0.24 to 0.6 by taking the arithmetic mean of their values at all grid points in the interval I = [0.24,0.6].For x > 0.6 (resp.x < 0.24), the boundary condition at x = 1 (resp.x = 0) influences the outcome, and, therefore, we leave values for x < 0.24 and x > 0.6 out of the averaging.The initial conditions for the CCEs Ẋ = c, (38) and (39) are set as follows: X(t 0 ) = X 0 is the location of the marginal tip density maximum, p(t 0 , x = X 0 , 0).We find X 0 = 0.18, set c(t 0 ) = c 0 = X 0 /t 0 .K(t 0 ) = K 0 is determined so that the maximum marginal tip density at t = t 0 coincides with the soliton peak.This yields K 0 = 196.
Using the soliton collective coordinates, we reconstruct the marginal vessel tip density and find its maximum value and the location thereof for all times t > t 0 .Figure 5 shows that the location of the soliton as predicted from the CCEs (38) and (39) compares very well with the marginal tip density obtained by ensemble averages (over 400 replicas) of the stochastic process during the 20 h time interval when boundaries do not affect soliton motion.There is a large discrepancy between the maximum marginal tip density as predicted by the soliton and by the stochastic process during the first 10 h of angiogenesis, which clearly marks the duration of the initial stage of soliton formation.After this stage, we note that the location of the maximum of the marginal tip density is very closely predicted by the location of the soliton peak as a function of time.However, the soliton density is noticeably wider than the ensemble averaged marginal density (see Figure 4).This is perhaps not so surprising.(i) The soliton is an approximate solution of Equation (28) for the marginal density that, in itself, is an approximation of (13) valid for large friction.A direct simulation of an overdamped Langevin model (1) with dv i (t) = 0 may provide better quantitative agreement with the soliton.Note that we do not have to select the velocity of the new tip branching from a given one in the overdamped case and that the injecting boundary condition (16) becomes −∂ p/∂x = 2βj 0 given by ( 19) at x = 0 [12].With a few exceptions [5][6][7]9], most models in the literature are overdamped.(ii) The soliton is an approximate solution for an unbounded slab and the influence of the boundaries has to be eliminated by selecting the interval I for spatial averaging.In view of Figure 5, we expect that the soliton approximation will improve if we increase the slab size so as to minimize the influence of boundaries.
In addition, and as it also happens in the simpler model of angiogenesis without haptotaxis [12], the maximum marginal tip density of different single realizations of the stochastic process (defined by (11) with N = 1) is well predicted by the position of the soliton peak density (see Figure 6).To explain this observation is an open problem.Comparison between the maximum value of p(t, x, 0) (calculated by ensemble average over 400 realizations) and its value as predicted by soliton collective coordinates.(a): of the maximum value of the marginal tip density (relative error smaller than 5.8% for times between 10 and 30 h); (b): evolution of the position of the maximum marginal tip density on [0, 1] (at t = 22 h, the absolute error is the space step in the numerical method, ∆x = 0.02; at t = 30 h, the error is 3∆x).
We have used an anastomosis coefficient Γ = 0.109.

Conclusions
In this paper, we have analyzed the influence of the haptotaxis mechanism in the first stages of tumor-induced angiogenesis following the Capasso-Morale model [7].Blood circulation, secondary branching, remodeling and pruning of capillaries [22,23] are not considered in this mesoscopic model, nor are microscopic features such as cell size, shape and cellular mechanics [4,18].The purpose of this paper is to show that the reduction of complex stochastic angiogenesis dynamics to the simpler soliton description of the angiogenic network allows incorporation of other transport mechanisms in addition to chemotaxis.As expected, including endothelial cell motion in an adverse fibronectin gradient through the extracellular matrix results in longer time scales for angiogenesis.As in previous work [9,10], there is a deterministic description associated to the full stochastic model: the density of active vessel tips satisfies an integrodifferential equation of Fokker-Planck type with a linear birth term and a nonlinear death (anastomosis) term.This equation is coupled to reaction-diffusion equations for the continuum fields (TAF, FN, MDE).Numerical simulations of the stochastic process show that, after an initial formation stage, the marginal tip density advances towards the tumor as a moving angiton whose longitudinal profile along the x-axis is soliton-like.Following the methodology explained in [11,12], we find an analytical formula for the soliton and collective coordinate equations that control shape and velocity of the soliton.When the soliton is about to reach the tumor, we need to study a different stage of soliton absorption by the tumor that is not contemplated in the present work.The detailed dynamics for TAF, FN and MDE as well as the chemotactic and haptotactic forces affecting vessel extension modify and influence the collective coordinate equations, which could be used to discuss possible anti-angiogenic strategies [28].The description of angiogenesis in terms of soliton dynamics is confirmed by comparison with numerical simulations of the full stochastic model.

Figure 1 .
Figure 1.Comparison between the vessel networks of two replicas after 36 h for the model with chemotaxis (a), and for the model with chemo and haptotaxis (b).The TAF level curves have also been depicted.

Figure 2 .
Figure 2. Effect of haptotaxis on the profile of p(t, x, 0) (calculated by ensemble average over 400 realizations) at times 24 h (a) and 36 h (b).Solid line: profile including haptotaxis, dashed line: profile without haptotaxis.

Figure 3 .
Figure 3. Normalized consumption of fibronectin given by f (t, x, y)m(t, x, y)/[ f (0, x, y)m(0, x, y)], cf.Equation (8), by the angiogenic network calculated for a single replica in panels (a) 24 h and (b) 36 h.Panels (c) and (d) correspond to (a) and (b) but now the consumption has been calculated by ensemble averages over 400 replicas.

Figure 4 .
Figure 4. Comparison between the profile of p(t, x, 0) (solid line, calculated by ensemble average over 400 realizations) and its reconstruction from the soliton collective coordinates (dashed line) at times: (a) 12 h; (b) 18 h; (c) 24 h; (d) 30 h.

Figure 5 .
Figure 5.Comparison between the maximum value of p(t, x, 0) (calculated by ensemble average over 400 realizations) and its value as predicted by soliton collective coordinates.(a): of the maximum value of the marginal tip density (relative error smaller than 5.8% for times between 10 and 30 h); (b): evolution of the position of the maximum marginal tip density on [0, 1] (at t = 22 h, the absolute error is the space step in the numerical method, ∆x = 0.02; at t = 30 h, the error is 3∆x).We have used an anastomosis coefficient Γ = 0.109.

Figure 6 .
Figure 6.Position of the soliton peak compared to that of the maximum marginal tip density for different replicas of the stochastic process.

Table 3 .
Dimensionless parameters in the hybrid model as given in