The statistical theory of the angiogenesis equations

Angiogenesis is a multiscale process by which a primary blood vessel issues secondary vessel sprouts that reach regions lacking oxygen. Angiogenesis can be a natural process of organ growth and development or a pathological induced by a cancerous tumor. A mean field approximation for a stochastic model of angiogenesis consists of partial differential equation (PDE) for the density of active tip vessels. Addition of Gaussian and jump noise terms to this equation produces a stochastic PDE that defines an infinite dimensional L\'evy process and is the basis of a statistical theory of angiogenesis. The associated functional equation has been solved and the invariant measure obtained. The results are compared to a direct numerical simulation of the stochastic model of angiogenesis and invariant measure multiplied by an exponentially decaying factor. The results of this theory are compared to direct numerical simulations of the underlying angiogenesis model. The invariant measure and the moments are functions of the Korteweg-de Vries soliton which approximates the deterministic density of active vessel tips.


Introduction
Angiogenesis is the process of cells organizing themselves into into blood vessels that grow from existing vessels and carry blood to organs and through tissue.It occurs in normal conditions of organ growth and regeneration, wound healing and tissue repair, and also in pathological conditions such as cancer, diabetes, rheumatoid arthritis or neovascular age-related macular degeneration.
Angiogenesis is driven by Vessel Endothelial Growth Factor (VEGF) and other pro-angiogenic proteins which are secreted by cells experiencing lack of oxygen (hypoxia).VEGF diffuses in the tissue, binds to extracellular matrix (ECM) components and forms a spatial concentration gradient in the direction of hypoxia.Once VEGF molecules reach an existing blood vessel, the latter walls open as a response and new vessel sprouts grow out of endothelial cells (ECs) off the vessel.Through the cellular signaling Notch process, VEGF activates the tip cell phenotype in ECs, which then grow filopodia with many VEGF receptors.The tip cells pull the other ECs (called stalk cells), open a pathway in the ECM, lead the new sprouts, and migrate in the direction of increasing VEGF concentration [26].Signaling and mechanical cues between neighboring ECs cause branching of new sprouts [28,31,48].Stalk cells in growing sprouts alter their shape to form a lumen (wall of the sprout) connected to the initial vessel that is capable of carrying blood [25].Sprouts meet and merge in a process called anastomosis to improve blood circulation inside the new vessels.Poorly perfused vessels may become thinner and their ECs, in a process that inverts angiogenesis, may retract to neighboring vessels leading to a more robust blood circulation [24].Thus, the vascular plexus remodels into a highly organized and hierarchical network of larger vessels ramifying into smaller ones [44].In normal processes of wound healing or organ growth, the cells inhibit the production of growth factors when the process is finished.In pathological angiogenesis, e.g., cancer, tumor cells lack oxygen and nutrients and produce VEGF that induces angiogenesis from a nearby primary blood vessel.The generated new vessel sprouts move and reach the tumor [23,15].Tumor cells continue secreting growth factors that attract more vessel sprouts and facilitate their expansion.
Together with experiments, many models spanning from the cellular to macroscopic scales try to understand angiogenesis; see the reviews [1,11,14,27,33,36,41,49].Early models consider reaction-diffusion equations for densities of cells and chemicals (growth factors, fibronectin, etc.) [32,13,18,19].They cannot treat the growth and evolution of individual blood vessels.Tip cell stochastic models of tumor induced angiogenesis are among the simplest ones for this com-plex multiscale process.Their basic assumptions are that (i) the cells of a blood sprout tip do not proliferate and move towards the tumor producing growth factor, and (ii) proliferating stalk cells build the sprout along the trajectory of the sprout tip.Thus tip cell models are based on the motion of single particles representing the tip cells and their trajectories constitute the advancing blood vessels network [42,43,1,36,33,7,27,11,45].Tip cell models describe angiogenesis over distances that are large compared with a cell size, thereby not incorporating descriptions of cellular and subcellular scales.These models are typically random: the motion of each tip typically includes directed Brownian motion, branching and anastomosis are birth and death processes, respectively.The random evolution of the tip cells may affect and be affected by reactiondiffusion equations for VEGF and other quantities (hybrid models).Other models contemplate motion of the tip cells on a lattice [1,36], and are related to cellular automata models [40,34,35].More complex models include tip and stalk cell dynamics, the motion of tip and stalk cells on the extracellular matrix outside blood vessels, shape and size changes of cells, signaling pathways and EC phenotype selection, blood circulation in newly formed vessels, and so on [2,30,46,41,3,47,27,31,39,49,4,11,48].The evolution of the blood vessels governed by these processes may be difficult to measure but various methods including laboratory experiments can be used to measure the statistical properties of an angiogenic network.These measurements of statistical quantities of the network can then be compared to simulations.
The authors have analyzed a system of stochastic differential equations plus birth and death processes, from [7,45,12], describing vessel growth and the associated PDE describing the evolution of the density of active vessel tips.In [9], they found that the growth of the tip of the blood vessel is well described by the KdV soliton.This leads to a control theory of angiogenesis that is currently being developed.In this paper, we will add the noise associated with branching and anastomosis, to the density equations, and develop the statistical theory of angiogenesis.Not surprisingly the solitons are found to play a major role in this theory and the noise in both branching and anastomosis is found to be multiplicative and depend only on the density.
The paper is organized in the following matter.In Section 2, we derive the density equations and formulated the noise in them.In Section 3, we explain the log-Poisson process arising from the noise in anastomosis and derive the geometric Brownian motion resulting from the noise in branching.In Section 4, we explain the statistical theory by the example of the stochastic Navier-Stokes Equation and derive the invariant measure of the stochastic angiogenesis equation.In Section 5, we compute the moments of the density and compare with simulations.In Section 6, we compute the structure function of the density.Section 7, contains a discussion.In Appendix A, we include a short summary of jumps and Lévy processes In Appendix B, we indicate how to calculate the moments of the density from numerical simulations of the angiogenesis stochastic process.

The Density Equation
We start with the equation for the density of active tip cells p in nondimensional form (1) from [7,45,12], where p(t, x) = p(t, x, v)dv, S(t, x) = t 0 p(s, x)ds, are the marginal density of p(t, x, v) and the density of stalk cells, respectively.δ v (v) is a Gaussian function having zero mean and small variance σ v .The density equation is derived from the Langevin equations for the blood vessel extension, see [7,45,12]: for t > T k , the random time when the kth tip, located at X k (t) and moving with velocity v k (t) appears as a consequence of another tip's branching.When an active tip arrives at a point that was occupied by another tip at a previous time, it disappears, which corresponds to anastomosis or loop formation [7,45,12].W k (t) denotes independent Brownian motion and the derivation uses Ito's lemma, where the Brownian noise √ σdW k (t) produces the Laplacian term σ 2 ∆ v p in the equation.The chemotactic force is where C is the VEGF concentration (called tumor angiogenic factor or TAF in tumor induced angiogenesis).In the hybrid stochastic model, the equation for the TAF density C(t, x) involves diffusion and consumption by the advancing tip cells [12], where N(t) is the number of active tips at time t and δ σ x are regularized delta functions: The source terms on the right hand side (RHS) of Eq. ( 1) arise from branching and anastomosis of active tips.The probability that a tip branches from one of the existing ones during an infinitesimal time interval (t,t + dt] is proportional to The velocity of the new tip that branches from tip i at time T i is selected out of a normal distribution, δ σ v (v − v 0 ), with mean v 0 and a narrow variance σ 2 v .The regularized delta function δ σ v (x) is given by Eq. ( 5) with σ x = σ y = σ v .
When using the deterministic description of Eq. ( 1), the C the tumor angiogenic factor satisfies the diffusive mean-field equation instead of Eq. ( 4).Here is the current density, see [12].Representative values of all involved dimensionless parameters can be found in Refs.[7,45,12].Density, marginal density and current density are the expected values of respectively, with respect to the stochastic process, as the regularizations of the delta functions disappear [45,12].Equations ( 1) and ( 7) are mean field approximations for these expected values and for the average TAF concentration.Moments of the density can be directly calculated from simulations of the stochastic process, as we shall see later.However, we want to build a theory of these moments by adding appropriate noise terms to the equation ( 1) and then analyzing the resulting stochastic PDE.The first term on the RHS of the equation (1) represents a (multiplicative) jump term and we add a multiplicative jump noise term p R h(z, v,t) N(dz, dt) associated with such jumps.The justification for this noise term is that tip branching and anastomosis are not deterministic processes, although they are represented in the density equation above as deterministic terms, in the mean field approximations.Thus we can expect noise to be associated to these terms and since both trip branching and anastomosis create jumps in the density, the correct form of the associated noise is multiplicative Poissonian jump noise.This gives the equation where α(C) = AC 1+C and b k t are i.i.d.Brownian motions.The first noise term is associated to velocity fluctuations through the Fourier coefficients < e ik•v , p >.Here h(z, v,t) = h b (z, v,t) + h a (z,t).We expect the noise to be small, ε << 1, in many cases but we will allow it to be as large as ε = 1, in the computations below.
denote the sizes of the jumps, associated with anastomosis and branching respectively, and N is the compensated number (of jumps) process, see [37].
The Fourier coefficients d k and c Thus the first and the last terms in the second line in the equation represent respectively the continuous (Brownian) and discrete (jump) noise.The Brownian terms are accompanied by a deterministic estimate for the large deviation (the d k s).We have made the Brownian terms as generic as possible by modeling noise in every (velocity) Fourier component of p.The stochastic partial differential equation ( 8) must be accompanied by vanishing boundary conditions on the plane x ∈ R 2 , and periodic boundary conditions in v ∈ T 3 .An initial condition for the density must also be specified.For the existence theory in the deterministic case, see [16] and [17], for convergence of positivity preserving numerical schemes, see [10].

The Log-Poisson Processes
An integrating factor can be found to simplify the density equation (8).It uses the Ito-Lévy formula and the geometric Lévy process.Let f (X k , v k ) = ln(q P ) and consider Ito-Lévy's formula (A.2) in Appendix A, where we have only taken the Poisson noise into account, β−1 , this means that neither h nor N t depend on z and the integral above reduces to a product.N t is the Poisson number process and − γ ln(b) β−1 its mean.The two parameters γ and b will be assigned values shortly.
see [6].Provided q P (0) = 1, solving for q gives This log-Poisson process is the integrating factor that we will use below.In [6], Example 1.5, it is shown how to compute the moments of a log-Poisson processes: 1. Suppose the Poisson process N k has the mean λ = − γ ln |k| β−1 (i.e., b = |k| above), k is the wave number, then it is straight-forward to compute the mean of the log-Poisson process |k| γ β N k .Namely, and we get the mean 2. Now we compute the nth moment E([|k| γ β N k ] n ) of the log-Poisson process |k| γ β N k above.By a similar computations as above, where λ is the mean from Part 1, therefore, Finally, we get that β−1 ) .(10) This will be the case for the example of the stochastic Navier-Stokes Equation below, we get a log-Poisson process for each wave number.

The Geometric Brownian motion gives the factor
3) for z in Appendix A, following Øksendal [38], where the initial conditions for the Brownian motions, starting at zero, are b k 0 .(Ito's formula produces the term (1/2)c k from the Brownian motion b k t that cancels the term −(1/2)c k in the corresponding large deviation.) is the drift coefficient, that we take to be negative.The initial k b k 0 , can either be positive or negative.

The Invariant Measure
The PDE (8) can be used to define an infinite dimensional Lévy process.The statistical properties such as the mean and moments of this process are determined by the invariant measure on function space associated with this process.We now review the theory of the Kolmogorov-Hopf equations determining the invariant measure.The first such equation that is a functional differential equation, where the derivatives are with respect to a function in a Banach space, was written down by Hopf [29], here φ is a bounded function of the Lévy process u, and Au is the deterministic part of the evolution equation determining u, •, • is the dual pairing in the Banach space.Hopf actually worked with the equation for the characteristic function, that is equivalent to the above equation, and he was looking for the invariant measure of the Navier-Stokes equation.Hopf's equation was reportedly solved by Kolmogorov that found that the invariant measure for the deterministic Navier-Stokes equation is disappointingly µ = δ(u), a delta function concentrated at the origin.The reason for this is that the Navier-Stokes equation is dissipative and all solutions eventually decay to the origin.
Da Prato and Zabczyk developed [21] a method that can be used to find the invariant measure for stochastic partial differential equation (SPDE) of the form where A is an operator on the Hilbert space on an n-torus L 2 (T n ).Here the c 1/2 j are n-vectors, their inner products converge where C is the trace class matrix having the c j s along the diagonal.The invariant measure producing the solution of this equation is Gaussian, where Q = ∞ 0 e tA Ce tA * dt is the variance, see [20].In [5] the invariant measure of the Navier-Stokes equation was computed using Stochastic Closure Theory, were the Navier-Stokes equation is split into a Reynolds Averages Navier-Stokes (RANS) equation and a stochastic Navier-Stokes (SNS) equation for the small scale flow.This is the equation where we have eliminated the pressure, using the incompressibility condition and the d k s terms are contributed by the large deviation principle, see [5].This equation can be written in the form The large deviation term in the noise will simply contribute a non-trivial mean to the infinite-dimensional Gaussian, but the last term that is multiplicative in u has to be treated differently.
In [5] the invariant measure of the stochastic Navier-Stokes equation was computed by applying Girsanov's theorem.This amounts to choosing the parameters γ = 2/3 = β and introducing an integrating factor based on the log-Poisson processes in the last section.This means that the Kolmogorov-Hopf equation becomes where C is the trace class matrix having the c j s along the diagonal, and D = [|k| 1/3 d k ] is an infinite vector.The invariant measure producing the solution of this equation is Gaussian, where Q = ∞ 0 e tA P t CP * t e tA * dt is the variance, and the mean is Tr(EI), where E = ∞ 0 e tA P t Ddt.D is the matrix with |k| 1/3 d k s along the diagonal and I is the matrix with Fourier components e k (x) on the diagonal.N k t in the integrating factor P t above, is the log-Poisson process and its law is: where m k t = − γ ln(|k|) β−1 is the mean of the Log-Poisson process for the kth wavenumber and is the probability of having exactly j jumps in the kth wavenumber.Thus we read the normal law with variance Q from the first and the last term of the Kolmogorov-Hopf equation, the mean from the second term and the Poisson law from the integrating factor.Then one can prove directly that ( 14) is the invariant measure of ( 12) and provides the solution of the Kolmogorov-Hopf equation ( 13), see [6].
We will now use the same reasoning to find the Kolmogorov-Hopf equation associated to (8) and solve it.First we write (8) as a stochastic partial differential equation, where pk =< p, e ik•v >, models velocity fluctuations.This shows that ( 16) is different from the stochastic Navier-Stokes equation above in that it contains two multiplicative noise terms and no additive noise.We write it in the form where we have set ε = 1.Now we find an integrating factor that is a product of two terms, a log-Poisson process as in the previous section, with the mean λ = m t = − γ ln(p t ) β−1 , and the geometric Brownian motion as in Equation (11).By Ito's formula Notice that the log-Poissonian is different from the stochastic Navier-Stokes equation where the jumps depend on the Fourier coefficients (k) of the solution, whereas for ( 16) the jumps depends on the whole solution p t , or are the same for all Fourier coefficients.Then the integrating factor becomes P t = q B q P and the Kolmogorov-Hopf equation reduces to Hopf's equation with an integrating factor The invariant measure of ( 16) is now a convolution of the Poisson distribution with a Gaussian, where the mean of the Gaussian is e t = ∑ k =0 (d k − 1 2 c k ) pk and the variance is The Gaussian is the distribution of the sum of the Fourier coefficients of p with respect to the velocity v, with mean e t and variance q t .The solution of the equation ( 17) can be written, see [6], as where m t = − γ ln(p t ) β−1 is the mean of the Log-Poisson process and P j = m j t e −m t j! is the probability of having exactly j jumps by the computation of E(q B ) and E(q P ) in Equations ( 11) and (10).We would like to take the limit t → ∞, to obtain the invariant measure where β−1 and e = e ∞ , q = q ∞ are the limits as t → ∞.However, this would give us the trival limit zero, because of the decay e −ndt due to the multiplicative Brownian noise.Thus we take the limit p t → p ∞ and consider the long time statistics to be determined by This can be thought of as an invariant measure multiplied by an exponentially decaying factor.It gives an excellent approximation to the real statistics for large times.

Comparison with Simulations and Theory
In [9], the authors showed that the density p evolves towards a 1D (Korteweg-de Vries) soliton profile for a slab geometry in which the primary vessel is the y axis and the tumor is centered at x = L, y = 0.A general 2D geometry is discussed in [12].For the 1D slab geometry, p is a product of a Maxwellian distribution in velocity, centered at v 0 , a transversal Gaussian function that approaches δ(y), and the soliton from [9], see [12].The justification of the Maxwellian distribution is that the source term in Equation ( 16) selects velocities in a small neighborhood of v 0 because they are the only velocities for which the birth term in the equation can compensate the anastomosis death term.Then a Chapman-Enskog expansion reduces Equation ( 16) to the equation in [9] that has the soliton solution, see [8].The factor δ(y) follows from a multiple scales method explained in [12].It is reasonable to expect the averaged density and its moments (19) to be functions of this soliton.This is in fact the case, but we now spell out the assumption that we make to establish this relation.The jump associated with Poisson process is: where 1 a and 1 b are one during branching and anastomosis respectively and zero otherwise.Now substituting in the soliton for p = a sech 2 (b(x − ct + x o )), we get that This expression ranges from 0 to 2(a/bc), with x and we choose the average S(∞) = (a/bc).Another way of determining this limit, is to let x + x o = ct and take the limit as x → ∞ and t → ∞.The (velocity) mean of the jump is then, the average of h over the velocity, This gives the average rate The averaged branching rate m b = 24.70 and the averaged anastomosis rate is m a = −22.49,see Table 1.However, the rates in Equation ( 19) depend on x (and t once we switch on the time evolution), so we must average them to get the average rates.We use the ergodic hypothesis to do this.Strictly speaking we do not have a non-trivial invariant measure to apply the ergodic theory, as discussed in the last section.However, we can ignore the decaying part q B of the solution and use (21) for the purpose of computing the averaged branching and anastomosis rates.We take the log of both sides of the equation from last section, and approximate the time averages of both sides by the averages with respect to (20) since E[log(q P )] ≤ log(E[q P ]) = 0, by Jensen's inequality, and E[q P ] = 1, by Equation (9).The first equality above is by the ergodic hypothesis, using the invariant measure (20), the second equality is by the mean value theorem, where p * is the value of p as some fixed velocity v * .We have set p = p ∞ q P and have ignored the decaying part q B .Finally, we approximate p * by the spatial mean of p A substitution of the soliton above into the integral gives . This allows to compute the exponents  25) to those computed from the stochastic simulations.Note that h a = − 1 2 ΓS(∞), while d 0 = 0. in Equation (19), using the values K = 266, c = 3, Γ = 0.32, µ = 8.5, used in the simulations, see [9], so B = 291.977and ln B = 5.677.For consistency, we have set α = 3.35.
Equation (19) now gives that assuming that the anastomosis and branching processes are independent, and using the values from Table 1, and since S(∞) = a/bc = 2KΓ + µ 2 = 15.57.The upshot of this is that unless we are in the traveling frame of the soliton all the moments are exceedingly small.However, in the traveling frame of the soliton ξ = x − ct + x 0 , the average is just the soliton itself whereas the second moment is and the nth moment is where ).
Now by the argument in the previous section, we get the moments In particular for the mean, we get a one parameter (d) family, whereas for the second moment, we get an expression depending on two parameters, d and ε, where instead of choosing the size of the perturbation in Equation ( 8) to be 1, we have chosen it to range between 0 and 1, as measured by ε, 0 ≤ ε ≤ 1.
Comparison of the theory and simulations is shown in Figure 1, for the second and the third moment at increasing times.It is reasonable to consider the decaying statistical quantities above instead of taking the limit as t → ∞.We are interested in times until a blood vessel (artery or vein) encounters another one or joins with a tumor.The decay coefficient is small d ≈ 1.8/hour so these vessels shrink slowly as they are elongated.The coefficient d 0 = 0, from the simulations.
There is a small discrepancy between the theoretical moments and the simulated ones at 24 hours, see Figure 1, both for the second and the third moment.This reason for this is well know.The soliton is not stable but there is a onedimensional subspace of translated solitons that is stable.This is called orbital stability, see [50].The moments are also orbitally stable so the theoretical moments are a small translate of the simulated ones, at 24 hours.This translation distance increases with time.

The Structure Functions
We compute the moments of the density differences p 1 − p 2 = p(ξ 1 ) − p(ξ 2 ) to probe the fine structure in angiogenesis, This means the we are looking at the difference of two solitons that differ only in their initial position.
The density differences δp We let ξ 1 = ξ 0 + l and ξ 2 = ξ 0 , where l is a lag variable.Then δp becomes by the same arguments as above, the second moment is for l small and the nth moment is for l small, where ).
These formulas then give the structure functions, as in the previous section.

Discussions
Considering the scaling exponents of the moments and the structure functions we see that the nth moment has a scaling n(1 − γ a − γ b ), with intermittency corrections γ a ( . This is different from the three-dimensional stochastic Navier-Stokes equation [5], where the moments are skewed Gaussians, but the structure functions have a scaling with intermittency corrections.The structure functions of ( 16) have the same scaling as the moments, however, moments of ( 16) are even functions and reach their maximum at ξ = 0, whereas the structure functions are odd functions of ξ 0 and have two maxima at ξ 0 = ±sech −1 ( 3/2).
We thus see that the statistical theory of ( 16) is very different from that of the stochastic Navier-Stokes equation.If we are in the traveling frame of the soliton, we see decaying soliton-like terms, given by the formulas above.Apart from these, the statistical quantities consist of small and rapidly decaying radiation terms.
We also see that a much simpler perturbation term, with continuous increments, of the density gives the same results, namely )dt − d 0 db t p + εp R h(z, v,t) N(dz, dt), (27) with b 0 = 1.Thus the (Brownian) noise, in vessel branching, only depends of the density (not all of its Fourier components) and is, in this aspect, similar to the noise in anastomosis, that also only depends on the jumps in the density.
Finally, the statistical theory of angiogenesis only applies to finite times, until a network of blood vessels is beneficially established in organs and recovering tissue, or malignantly connected a cancerous tumor.
Summarizing, we have proposed a statistical theory of angiogenesis by adding appropriate noise terms to the equation for the density of active tip cells.The shape of the noise terms mimic the birth and death processes of tip branching and anastomosis, respectively, plus some other terms coming from Brownian motion.Thus, our start point is a stochastic PDE with Gaussian and Poisson noises that define an infinite dimensional Lévy process.We have solved the associated functional equation by finding appropriate integrating factors and therefore obtained the invariant measure.The result is an invariant measure multiplied by an exponentially decaying factor.By comparing theory and numerical simulations of the underlying angiogenesis model, we have obtained the appropriate values of the parameters involved in the stochastic PDE.We have found that the invariant measure and the moments are functions of the soliton which approximates the deterministic density [9,8].
It is illustrative to use the (A.2) to solve a differential equation.We now do so creating the geometric Lévy process: We solve the differential equation with respect to t and exponentiating, we get that 2 )t+αb t + t 0 R ln(1+h(s,z)) N(dz,ds)+ t 0 R (ln(1+h(s,z))−h(s,z))m(dz)ds} .(A.3) This process is called the geometric Lévy process.The first two terms in the exponent correspond to the Ito process and the last two terms (integrals) are contributed by the jump (Lévy) process.
and the b j t are independent Brownian motion accompanying each Fourier coefficient.The corresponding Kolmogorov-Hopf equation is

Figure 1 :
Figure 1: A comparison of theory and simulations for the second (first three) and the third (last three) moments of the stochastic angiogenesis equation (16).The times are 16 hours (left), 20 hours (middle) and 24 hours (right), from left to right.