Abstract
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 one induced by a cancerous tumor. A mean-field approximation for a stochastic model of angiogenesis consists of a partial differential equation (PDE) for the density of active vessel tips. Addition of Gaussian and jump noise terms to this equation produces a stochastic PDE that defines an infinite-dimensional Lévy 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 of this theory are compared to direct numerical simulations of the underlying angiogenesis model. The invariant measure and the moments are functions of a Korteweg–de Vries-like soliton, which approximates the deterministic density of active vessel tips.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
Angiogenesis is the process of cells organizing themselves 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 the Vessel Endothelial Growth Factor (VEGF) and other pro-angiogenic proteins that 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 notch signaling 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 (Gerhardt et al. 2003). Signaling and mechanical cues between neighboring ECs cause branching of new sprouts (Hellström et al. 2007; Jolly et al. 2015; Vega et al. 2020). 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 (Gebala et al. 2016). 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 (Franco et al. 2015). Thus, the vascular plexus remodels into a highly organized and hierarchical network of larger vessels ramifying into smaller ones (Szymborska and Gerhardt 2018). 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 (Folkman 1971, 1974; Carmeliet 2005). 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 (Anderson and Chaplain 1998; Bonilla et al. 2019; Byrne 2010; Heck et al. 2015; Mantzaris et al. 2004; Plank and Sleeman 2004; Scianna et al. 2013; Vilanova et al. 2017). Early models consider reaction–diffusion equations for densities of cells and chemicals (growth factors, fibronectin, etc.) (Liotta et al. 1977; Byrne and Chaplain 1995; Chaplain and Stuart 1993; Chaplain 1995). 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 complex multiscale process. Their basic assumptions are that (i) the cells of a blood sprout tip do not proliferate and move toward 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 vessel network (Stokes and Lauffenburger 1991; Stokes et al. 1991; Anderson and Chaplain 1998; Plank and Sleeman 2004; Mantzaris et al. 2004; Bonilla et al. 2014; Heck et al. 2015; Bonilla et al. 2019; Terragni et al. 2016). 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 reaction–diffusion equations for the VEGF and other quantities (hybrid models). Other models contemplate motion of the tip cells on a lattice (Anderson and Chaplain 1998; Plank and Sleeman 2004) and are related to cellular automata models (Pillay et al. 2017; Martinson et al. 2020, 2021). 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 (Bauer et al. 2007; Jackson and Zheng 2010; Travasso et al. 2011; Scianna et al. 2013; Bentley et al. 2014; Van Oers et al. 2014; Heck et al. 2015; Jolly et al. 2015; Perfahl et al. 2017; Vilanova et al. 2017; Bernabeu et al. 2018; Bonilla et al. 2019; Vega et al. 2020). 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 and experiments.
The authors have analyzed a system of stochastic differential equations plus birth and death processes (Bonilla et al. 2014; Terragni et al. 2016; Bonilla et al. 2020), describing vessel growth, and the associated PDE describing the evolution of the density of active vessel tips. In Bonilla et al. (2016b), they found that the growth of the tips of the blood vessels is well described by a Korteweg–de Vries-like (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 equation and develop the statistical theory of angiogenesis. Not surprisingly the soliton is found to play a major role in this theory, and the noise in both branching and anastomosis is found to be multiplicative and to depend only on the density of active tip cells.
The statistical theory is a crucial tool to determine the properties of the angiogenic network of veins. These cannot be found by simulating or measuring a single vein and to determine all the properties of the network we need the moments of the density of active tip cells. It may also be necessary, for the finer structure of the network, to compute the structure functions that are the moments of the difference between densities separated by a spatial lag variable. In this paper, we develop the tools to compute the moments and structure functions of the density.
The paper is organized in the following matter. In Sect. 2, we outline the derivation of the density equation for the active blood vessel tips, from previous papers, and formulate the macroscopic noise in them. In Sect. 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. These are used in the next section to compute the moments of the density. In Sect. 4, we explain the statistical theory and derive the invariant measure of the stochastic angiogenesis equation. The moments of the density are computed using the invariant measure and the results from Sect. 4. In Sect. 5, we explain the geometry and the form of the soliton and compute the coefficients that the moments of the density depend on. Then the computed moments are compared with simulations. In Sect. 6 we compute the structure function of the density. Section 7 contains a discussion. In each section, the results are frequently compared to the analogous results for the stochastic Navier–Stokes equation, the only nonlinear stochastic PDE that this statistical theory has been previously developed for, see Birnir (2013a, 2013b). In “Appendix A,” we include a short summary of jump 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. In “Appendix C,” we present a slight extension of the mean-value theorem that is used in Sect. 5.
2 The Density Equation
In this section, we outline the derivation of Eq.�(1) describing the evolution of the density for the active tips of blood vessels, see Bonilla et�al. (2014), Terragni et�al. (2016), Bonilla et�al. (2020). Then, we formulate the macroscopic noise associated with tip branching and anastomosis and add it to the density equation. This gives a nonlinear stochastic PDE, namely Eq.�(9) below, for the noisy evolution of the tip density.
There is noise both in the stochastic equations on the microscopic level and in the density equation on the macroscopic level. The origin of the noise in both equations is explained below. It is natural to assume that the noise, with continuous increments, in the density equations depends on the velocity components, only after the computations does it become clear that it only depends on the density.
We start with the equation for the density of active tip cells, \(p(t,x,v) \in L^2({\mathbb {R}}^+;\ {\mathbb {R}}^2,\ {\mathbb {R}}^2)\), where p is a function of time t, location x and velocity v. In nondimensional form, it is
obtained from (Bonilla et�al. 2014; Terragni et�al. 2016; Bonilla et�al. 2020), where
are the marginal density of p(t, x, v) and the density of stalk cells, respectively, and \(\Gamma \) is a constant. Here, \(\delta _v(v)\) is a delta function, regularized with a Gaussian
with zero mean and small standard deviation \(\sigma _v\).Footnote 1 The density equation is derived from the Langevin equations for the blood vessel extension, see Bonilla et al. (2014), Terragni et al. (2016), Bonilla et al. (2020):
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. The tip position \(X^k\) is slave to the velocity that is sensitive to (microscopic) noise. This is why the Brownian noise, with variation \(\sigma \), appears in the \(v^k\) equation. 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 (Bonilla et al. 2014; Terragni et al. 2016; Bonilla et al. 2020). In Eq. (3), \(W^k(t)\) denotes independent Brownian motions. The derivation of (1) uses Ito’s lemma, where the Brownian noise \(\sqrt{\sigma } dW^k(t)\) produces the Laplacian term \(\frac{\sigma }{2}\Delta _v p\) in the density equation. The chemotactic force is
where C is the VEGF concentration (called tumor angiogenic factor or TAF in tumor-induced angiogenesis), while \(\Gamma _1\), \(\delta _1\) and q are constants. In the hybrid stochastic model, the equation for the TAF density \(C(t,x)\in C^0({\mathbb {R}}^+;C^2({\mathbb {R}}^2))\) (functions continuous in t and twice continuously differentiable in x) involves diffusion and consumption by the advancing tip cells (Bonilla et�al. 2020),
where N(t) is the number of active tips at time t and \(\delta _{\sigma _x}\) is a regularized delta function
The source terms in 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 \(\sum _{i=1}^{N(t)}\alpha (C(t,X^i(t))){d}t\), with
where A is a positive constant. The velocity of the new tip that branches from the ith tip at time \(T^i\) is selected out of a normal distribution, \(\delta _{\sigma _v}(v-v_0)\), with mean \(v_0\) and a narrow variance \(\sigma _v^2\). The regularized delta function \(\delta _{\sigma _\textbf{v}}(\textbf{v})\) is given by Eq.�(6) with \(\sigma _x=\sigma _y=\sigma _v\).
When using the deterministic description of Eq.�(1), the TAF satisfies the diffusive mean-field equation
with \(C(t,x)\in C^0({\mathbb {R}}^+;C^\infty ({\mathbb {R}}^2))\), instead of Eq.�(5). Here
is the current density, see Bonilla et�al. (2020). Representative values of all involved dimensionless parameters can be found in Refs. Bonilla et�al. (2014, 2020), Terragni et�al. (2016).
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 (Terragni et�al. 2016; Bonilla et�al. 2020). Equations (1) and (8) 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 (macroscopic) noise terms to Eq.�(1) and then analyzing the resulting stochastic PDE. The first term in the RHS of Eq.�(1) represents a (multiplicative) jump term; hence, we add a multiplicative jump noise term
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 approximation. Thus we can expect noise to be associated with these terms, and since both tip branching and anastomosis create jumps in the density, the correct form of the associated noise is a multiplicative Poissonian jump noise. Notice that the noise exists on many levels and the terms that we are adding to Eq.�(1) represent macroscopic noise, which is distinct from the microscopic noise in Eq.�(3). This gives the equation
where \(\alpha (C)\) is given by Eq.�(7) and \(b_t^k\) are i.i.d. Brownian motions. The first noise term is associated with velocity fluctuations, through the (velocity) Fourier coefficients \(<e^{ik\cdot v},p>=\int _{{\mathbb {R}}^2}\int _{{\mathbb {R}}^2}e^{ik\cdot v}pdvdx\), of the spatial mean. Here the cumulative jump size is \(h(z,v,t)=h_b(z,v,t)+h_a(z,t)\). We expect the noise to be small, \(\varepsilon<< 1\), in many cases, but we will allow it to be as large as \(\varepsilon =1\), in the computations below. In addition,
and
denote the sizes of the jumps, associated with anastomosis and branching, respectively, and \({{\bar{N}}}\) is the compensated number (of jumps) process, see �ksendal and Sulem (2005). The (velocity) Fourier coefficients \(d_k\) and \(c_k^{1/2}\) are arbitrary, but they are summable, \(\sum _{k\ne 0} |d_k| < \infty \), \(\sum _{k\ne 0} |c_k^{1/2}| < \infty \). The \(c_k^{1/2}\) are also square summable, \(\sum _{k\ne 0} c_k < \infty \). Thus, the first and last terms in the second line in Eq.�(9) represent the continuous (Brownian) and the discrete (jump) noise, respectively. 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. This is necessary because the velocity of a tip may be sensitive to noise in infinitely many directions in the Hilbert space of the velocity. Stochastic PDE (9) must be accompanied by vanishing boundary conditions on the plane \(x \in {\mathbb {R}}^2\) and vanishing boundary conditions on the velocity plane \(v \in \mathbb {R}^2\). This makes the Hilbert space be \(L^2({\mathbb {R}}^+;\ \Omega ,\ {\mathbb {T}}^3)\), where \(\Omega \subset {\mathbb {R}}^2\) is a rectangle with Dirichlet boundary conditions, see Bonilla et al. (2016b). In the remainder of the paper, we will analyze (weak or strong) solutions of Eq. (9) on this space, with an appropriate initial density. For the existence theory in the deterministic case, see Carpio and Duro (2016) and Carpio et al. (2017), for convergence of positivity preserving numerical schemes, see Bonilla et al. (2018).
3 The Log-Poisson Process and Geometric Brownian Motion
The noise terms in Eq. (9) give rise to log-Poisson processes and geometric Brownian motions. In this section we show how to compute their moments. This is then used in the subsequent section to solve the (Kolmogorov–Hopf) equation for the invariant measure and compute the moments of the density p.
An integrating factor can be found to simplify density Eq. (9). It uses the Ito–Lévy formula and the geometric Lévy process. Let \(g(X^k,v^k) = \ln (q_P)\) and consider Ito–Lévy’s formula (A.2) in “Appendix A,”
Here, if N(z, t) is the number process of the Lévy process \(\eta (t)\), then \(m(U) =\int _U E(N(\hbox {d}z,1))\) is the Lévy measure of \(\eta (t)\), and \({{\bar{N}}}(\hbox {d}z,\hbox {d}t) = N(\hbox {d}z,\hbox {d}t)-m(\hbox {d}z)\hbox {d}t\) if \(|z|<R\) or \({{\bar{N}}}(\hbox {d}z,\hbox {d}t) = N(\hbox {d}z,\hbox {d}t)\) if \(|z|\ge R\) is the compensated jump measure of \(\eta (t)\); see Birnir (2013b), Section 1.5. Now, suppose that we only take the Poisson noise into account,
Then, if \(h=\beta -1\) and the mean of the Poisson number process \(N_t=N(z,t)\) is \(E(N_t)=-\frac{\gamma \ln (b)}{\beta -1}\), this means that neither h nor \(N_t\) depends on z and the integral above reduces to a product. The two parameters \(\gamma \) and b will be assigned values shortly. Thus, the previous formula becomes
see Birnir (2013b). Provided \(q_P(0)=1\), solving for \(q_P\) gives
This log-Poisson process is the integrating factor that we will use below. In Birnir (2013b), Example 1.5, it is shown how to compute the moments of a log-Poisson process.
-
1.
Suppose that the Poisson process \(N_k\) has the mean \(\lambda = -\frac{\gamma \ln |k|}{\beta -1}\) (i.e., \(b=|k|\) above), where k is the wave number. Then it is straightforward to compute the mean of the log-Poisson process \(|k|^\gamma \beta ^{N_k}\):
$$\begin{aligned} E(|k|^\gamma \beta ^{N_k})=\sum _{j=0}^\infty |k|^\gamma \beta ^j \frac{\lambda ^j}{j!} e^{-\lambda }=|k|^\gamma \sum _{j=0}^\infty \frac{(\beta \lambda )^j}{j!} e^{-\lambda }=|k|^\gamma e^{(\beta -1)\lambda }. \end{aligned}$$Thus,
$$\begin{aligned} \ln [E(|k|^\gamma \beta ^{N_k})]= \gamma \ln |k|+(\beta -1)\lambda =\gamma \ln |k|-(\beta -1)\frac{\gamma }{\beta -1}\ln |k| =0, \end{aligned}$$and we get the mean
$$\begin{aligned} E(|k|^\gamma \beta ^{N_k})=1. \end{aligned}$$(10) -
2.
Now, we compute the nth moment \(E([|k|^\gamma \beta ^{N_k}]^{n})\) of the log-Poisson process \(|k|^\gamma \beta ^{N_k}\). By a similar computation as above,
$$\begin{aligned} E([|k|^\gamma \beta ^{N_k}]^{n})=|k|^{n\gamma } e^{(\beta ^{n}-1)\lambda }, \end{aligned}$$where \(\lambda \) is the mean from Part 1. Therefore
$$\begin{aligned} \ln [E([|k|^\gamma \beta ^{N_k}]^{n})]=n\gamma \ln |k|+(\beta ^{n}-1)\lambda =\!\left( n-\frac{\beta ^{n}-1}{\beta -1}\right) \!\gamma \ln |k|. \end{aligned}$$Finally, we get that
$$\begin{aligned} E([|k|^\gamma \beta ^{N_k}]^{n})=|k|^{\gamma \left( n-\frac{\beta ^{n}-1}{\beta -1}\right) }. \end{aligned}$$(11)This computation is similar for the stochastic Navier–Stokes equation, the difference being that there we get a log-Poisson process for each wave number, see Birnir (2013a).
-
3.
The geometric Brownian motion gives the factor
$$\begin{aligned} q_B=\exp \!\left\{ \sum _{k\ne 0}\left[ \left( d_k-\frac{1}{2}c_k\right) t+c_k^{1/2}b_t^k\right] \right\} \!, \end{aligned}$$see formula (A.3) for \(Z_t\) in “Appendix A,” following Øksendal (Øksendal 2003). Then,
$$\begin{aligned} E(q_B)=\exp \!\left\{ \sum _{k\ne 0} d_k\ t\ + \sum _{k\ne 0}c_k^{1/2}b_0^k \right\} \!=\exp {(-d\ t-d_0)}, \end{aligned}$$(12)where the initial conditions for the Brownian motions, starting at zero, are \(b_0^k\) (note that Ito’s formula produces the term \((1/2)c_k\) from the Brownian motion \(b_t^k\) that cancels the term \(-(1/2)c_k\) in the corresponding large deviation). Thus, \(-d = \sum _{k\ne 0} d_k\) is the drift coefficient that we take to be negative. The initial condition \(d_0=-\sum _{k\ne 0}c_k^{1/2}b_0^k\) can be either positive or negative.
4 The Invariant Measure
In this section, we give a brief outline of the theory of the Kolmogorov–Hopf equation determining the invariant measure of stochastic PDEs. We first discuss its derivation for the Navier–Stokes equation where it originated. Then, we derive the Kolmogorov–Hopf equation for stochastic density Eq. (9) and solve it. This determines the invariant measure in Eq. (16). We end the section by computing the moments \(\langle p^n \rangle \) of the density (cf Eq. (20)), where the average uses the invariant measure.
Stochastic PDE (9) 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 the 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 (1952),
Here, \(\phi \) is a bounded function of u and A is a partial differential operator. The deterministic evolution equation determining u is
and \(\langle \cdot , \cdot \rangle \) is the dual pairing in the Banach space. Hopf actually worked with the equation for the characteristic function, which 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 who found that the invariant measure for the deterministic Navier–Stokes equation is disappointingly \( \mu = \delta (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 a method (Da Prato and Zabczyk 1996) that can be used to find the invariant measure for a stochastic PDE of the form
where A is an operator on the Hilbert space on an n-torus \(L^2({\mathbb {T}}^n)\). Here, the \(c_j^{1/2}\) are n-vectors, their inner products converge, \(\sum _{j \in {\mathbb {Z}}^n} c_j \le \infty \), with \(c_j =c_j^{1/2}\cdot c_j^{1/2}\), and the \(b^j_t\) are independent Brownian motions accompanying each Fourier coefficient. The corresponding Kolmogorov–Hopf equation is
where C is the trace class matrix having the \(c_j\)s along the diagonal. This is a functional differential equation: \(\phi \) is a bounded function of the solution u of the PDE, and the gradient \(\nabla _u\) and the Laplacian \(\Delta _u\) are with respect to the function u. The invariant measure producing the solution of the Kolmogorov–Hopf equation is an infinite-dimensional Gaussian,
where \(Q= \int _0^\infty e^{tA}Ce^{tA^*}\hbox {d}t\) is the variance, see Da Prato (2006).
We will now find the Kolmogorov–Hopf equation associated with Eq. (9) and solve it. First we write Eq. (9) as a stochastic PDE, namely
where \({\hat{p}}_k =\langle e^{ik\cdot v},p\rangle \) models velocity fluctuations and S is evaluated as the mean of p, assuming ergodicity, as in Eq. (23) below. This shows that Eq. (13) is different from the stochastic Navier–Stokes equation, see Birnir (2013a), since it contains two multiplicative noise terms and no additive noise. We write it in the form
where A is a partial differential operator and we have set \(\varepsilon =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 \(\lambda =m_t = -\frac{\gamma \ln (p_t)}{\beta -1}\), and the geometric Brownian motion
By Ito’s formula (cf Eq. (12)) we get
Notice that the log-Poissonian is different from the stochastic Navier–Stokes equation, where the jumps depend on the Fourier coefficients of the solution, see Birnir (2013a), whereas for Eq. (13) the jumps depend on the whole solution \(p_t\) or are the same for all Fourier coefficients. Then, the integrating factor becomes
and the Kolmogorov–Hopf equation reduces to Hopf’s equation with an integrating factor,
The invariant measure of Eq.�(13) is now a convolution of the Poisson distribution with an infinite-dimensional Gaussian, namely
where the mean of the Gaussian is \(e_t=\sum _{k\ne 0} (d_k-\frac{1}{2}c_k){{\hat{p}}_k}\) and the variance is \(q_t =\sum _{k\ne 0}c_k{{\hat{p}}}_k^2\). Indeed, the Gaussian distribution depends on the sum of the (velocity) Fourier coefficients of p. On the other hand, \({\mathbb {P}}_{N_\infty }(\cdot )\) is the limit as \(t\rightarrow \infty \) of the log-Poisson law.
The solution of Eq. (15) can be written see Birnir (2013b)) as a Markov semigroup \(R_t\) acting on a bounded function \(\phi \), on the Hilbert space \(\mathcal {H}\) [for example a Sobolev space based on \(L^2\) as in Birnir (2013b)] containing p,
where \(m_t = -\frac{\gamma \ln (p_t)}{\beta -1}\) is the mean of the log-Poisson process, \({\mathcal {P}}_j = \frac{m_t^j e^{-m_t}}{j!} \) is the probability of having exactly j jumps, \(N^j_\infty = N^j\). Note that \(R_t\) leaves the invariant measure in Eq. (19) below invariant. For \(\phi = p^n\), we get
by the computation of \(E(q_B)\) and \(E(q_P)\) in Eq. (12) and (11), where \(e^{At}\) denotes the solution operator generated by the partial differential operator A and \({{\bar{p}}}_t\) is the averaged solution of Eq. (14). We would like to take the limit as \(t \rightarrow \infty \), to obtain the invariant measure
where \(m_\infty = -\frac{\gamma \ln (p_\infty )}{\beta -1}\) and \(e= e_\infty ,\, q=q_\infty \) are the limits as \(t \rightarrow \infty \). However, this would give us the trivial limit zero, because of the decay \(e^{-ndt}\) due to the multiplicative Brownian noise. In Bonilla et al. (2016b), it was shown that the density evolves toward a soliton-like solution as t becomes large. Thus, we take the limit \(p_\infty =p(\xi )\), with \(\xi =x-ct\), to be the soliton in its traveling frame and consider the long-time statistics to be determined by
This can be thought of as the moments of an invariant measure multiplied by an exponentially decaying factor. It gives an excellent approximation to the real statistics for large times. In other words, the infinite time limit of \(p_t\) is trivial as explained above, but instead \(p_\infty \) refers the soliton \(p(\xi )\) in its traveling frame of reference. It is the correct long-time average.
5 Comparison between Simulations and Theory
In this section, we will compare the moments of the density computed as in Eq. (20) with the moments simulated from the stochastic model of angiogenesis. But first we discuss the form of the soliton that approximates the density of the active vessel tips and then we explain how the parameters in Eq. (20) are computed.
In Bonilla et al. (2016b), the authors showed that the density p in Eq. (1) evolves toward a 1D (Korteweg–de Vries-like) soliton profile for a slab geometry in which a primary vessel is located at \(x = 0\) and a tumor is centered at \(x=L\), \(y=0\), see Fig. 1. A one-dimensional slab is the x-axis from 0 to L. A two-dimensional slab in the \(x-y\) plane is the rectangle with length L and width 3L, see Fig. 1. A general 2D geometry is discussed in Bonilla et al. (2020). For the 2D slab geometry, p is a product of a Maxwellian distribution in velocity, centered at \(v_0\), a transversal Gaussian function that approaches \(\delta (y)\) and the soliton from Bonilla et al. (2016b), namely
See Bonilla et al. (2020). The soliton given by Eq. (22) is a 1D traveling wave solution of Eq. (1) in which only dominant terms (time derivative of the density, constant convection and source terms) are kept (Bonilla et al. 2016b, 2020), see Fig. 1. Here K and c are slowly varying functions of t that fix the size and velocity of the soliton. They are called collective coordinates, depend on all nondominant terms in Eq. (1) and satisfy ordinary differential equations (Bonilla et al. 2016a). All other parameters in Eq. (22) are related to the underlying angiogenesis model. The justification of the Maxwellian distribution is that the source term in Eq. (13) 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 Eq. (13) to the equation in Bonilla et al. (2016b) that has the soliton solution, see Bonilla et al. (2016a). This reduction is a side issue for our discussion and is not repeated here. It is explained in great detail in Bonilla et al. (2016a). The factor \(\delta (y)\) follows from a multiple scales method explained in Bonilla et al. (2020).
It is reasonable to expect the averaged density and its moments in Eq. (18) to be functions of the soliton in Eq. (22), since it is the only structure that survives in the numerical simulations for large times. This is in fact the case, but we now spell out the assumption that we need to establish this relation. The jump associated with the Poisson process is
where \(\mathbbm {1}_b\) and \(\mathbbm {1}_a\) are one during branching and anastomosis, respectively, and zero otherwise. Now, substituting for the soliton \(\tilde{p} = a\,\text{ sech}^2(b(x-ct+x_0))\), we get that
The range of this expression is (0, 2(a/bc)) for positive values of t, and we have chosen the average \(S(\infty )=a/bc\). Another way of determining this value is to let \(x+x_0=ct\) and take the limit as \(x \rightarrow \infty \) and \(t \rightarrow \infty \). The mean of the jump is the integral of h over the velocity,
This gives the averaged rate
where \(\overline{\alpha }\) is an equilibrium concentration, see the derivation below. Now, the average denoted by a bar indicates the average over the invariant measure from last section. This will be equal to the ensemble average discussed in “Appendix B.”
The averaged branching rate \({\overline{m}}_b\) and the averaged anastomosis rate \({\overline{m}}_a\) are reported in Table 1. Similarly, \({\overline{h}}_a\) and \({\overline{h}}_b\) denote the averaged jumps in anastomosis and branching, respectively, measured in simulations (see Table 1). However, the rates in Eq. (18) depend on x (and t once we switch on the time evolution), so we must average them to get the averaged rates. We perform this averaging. Strictly speaking we do not have a trivial invariant measure to apply the ergodic theory, as discussed in the last section. However, we can ignore the decaying part \(e^{-dt-d_0}\) of the solution and use the \(p_\infty \) part of Eq. (20) for the purpose of computing the averaged branching and anastomosis rates. Recalling that the mean of the log-Poisson process in Eq. (17) is
where \(h = \beta -1\) is the jump, we use the averaged jumps for both the branching and anastomosis processes. Then, if we can evaluate the average of the mean with respect to the invariant measure, we get an equation for the parameter \(\gamma \) in both cases. This is done by averaging the mean with respect to the infinite-dimensional Gaussian measure. Indeed, it is already the mean with respect to the law of the log-Poisson process, so the Poissonian part of the invariant measure can be omitted. This gives
by the slight extension of the mean value theorem in “Appendix C,” where \(p^*\) is a finite positive value. We do not know what the value \(p^*\) is, but a reasonable guess is to approximate \(p^*\) by the spatial mean of the soliton \({\tilde{p}}\) in Eq. (22), namely
A substitution of the soliton in Eq. (22) into the integral gives
where \(B= \frac{2(\sqrt{2K\Gamma +M^2})c}{\Gamma }\). This allows to compute the exponents
in Eq. (18).
The values of the parameters \({{\overline{m}}}_b\), \({{\overline{m}}}_a\), \({{\overline{h}}}_b\) and \({{\overline{h}}}_a\) are given in Table 1 for the second and third moments, and the considered time instants. Note that the values of \({{\overline{m}}}_b\) and \({{\overline{m}}}_a\) can be estimated directly from the stochastic simulations, those of \({{\overline{h}}}_a\) come from \(\overline{h}_a = -\Gamma S(\infty ) = -\Gamma a/bc = -\sqrt{2K\Gamma +M^2}\), while the values of \({{\overline{h}}}_b = \overline{\alpha }\) [and d, also reported in Table 1, cf Eq. (18)] are selected as to fit the theoretical moments in Eq. (26) to the numerical simulations of the angiogenesis model (see “Appendix B”). On the other hand, values of the collective coordinates are computed, for each time instant, as detailed in Bonilla et al. (2016a). Finally, the other parameter values are set to \(\Gamma = 0.09\), \(M = 7.78\) and \(F_x = 0.08\) (see Bonilla et al. 2016a). It is important to remark that, while the numerical values of \(\overline{h}_a\) are negative since the anastomosis is a death process, a correct fitting is obtained by considering the actual size of the anastomosis jumps, namely \(|\overline{h}_a|\).
Equation (18) now gives that
for t sufficiently large, assuming that the anastomosis and branching processes are independent, where \(p_\infty (\xi )\) is simply the soliton in Eq. (22). Here
and
The upshot of this is that, unless we are in the traveling frame of the soliton, all the moments are exceedingly small. The reason is that the soliton decays exponentially in space, so all quantities are going to be small away from the center of the soliton. Thus the soliton decays exponentially as \(\xi \rightarrow \pm \infty \) and has significant values only around \(\xi = 0\). However, in the traveling frame of the soliton, \(\xi = x-ct+x_0\), the average is just the soliton itself,
and the nth moment is
Here
Now, by the argument in the previous section, we get the moments
where the angle bracket denotes the average with respect to the invariant measure in Eq. (19). In particular, for the mean we get a one parameter (d) family,
Comparison between theory and simulations, discussed in “Appendix B,” is shown in Fig. 2, for the second and third moments. Results are given in the spatial domain \(0 \le x \le 1\) at increasing time instants.
It is reasonable to keep the exponentially decaying factor above instead of taking the limit as \(t \rightarrow \infty \). Indeed, 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 (see Table 1), so these vessels shrink slowly as they are elongated. The coefficient \(d_0 = 0\) is set for all simulations.
There is a small discrepancy between the theoretical moments and the simulated ones at 24 h (see Fig. 2). The reason for this is well-known. The soliton is not stable, but there is a one-dimensional subspace of translated solitons that is stable. This is called orbital stability, see Weinstein (1987). The moments are also orbitally stable, so the theoretical moments are a small translation of the simulated ones, at 24 h. The translation distance increases with time.
6 The Structure Functions
The moments in Eq. (20) are called the 1-point statistics of the statistical theory (of angiogeneses). The 2-point statistics, given by the structure function below, sometimes allow one to probe a finer structure. This is the case for example in turbulence (Birnir 2013a). We will compute the structures functions below to see if this also works for angiogenesis.
We compute the moments of the density difference \(p_1-p_2=p(\xi _1)-p(\xi _2)\) to probe the fine structure in angiogenesis, with \(\xi _i = x-ct+x_i^0,\ i =1,2\). These moments are called the structure functions. This means that we are looking at the difference of two solitons that differ only in their initial position.
The density difference \(\delta p = p_1-p_2\) is
We let \(\xi _1=\xi _0+l\) and \(\xi _2=\xi _0\), where l is a lag variable. Then \(\delta p\) becomes
By the same arguments as above, the nth moment is
for l small, where
These formulas then give the structure functions
as in the previous section, where the angle bracket denotes the average with respect to the invariant measure in Eq. (19). It equals the ensemble average discussed in “Appendix B.”
7 Discussions
We have developed a statistical theory for stochastic angiogenesis equation (9) and compared the moments of the density, solving (9), to simulations of the moments. The moments of the density turn out to be functions of the soliton, that approximates the deterministic density (Bonilla et al. 2016a, b), that is a solution of deterministic Eq. (1). The mean of the stochastic density is computed to be the soliton itself, see Eq. (25). However, for the higher \(\ge 2\) moments, the jumps in branching and anastomosis also come into play. These jumps, that are computed from simulations in Sect. 5, determine the scaling \(\gamma _b\) and \(\gamma _a\) and intermittency \(\beta _b\) and \(\beta _a\) parameters.
All of these moments and the structure functions are functions of the soliton but multiplied by a slowly decaying factor. This reflects the numerical observations: The density approaches a slowly decaying soliton for large times. It turns out to be possible to separate the slow decay from an invariant motion and compute the invariant measure determining this motion.
Considering the scaling exponents of the moments
We see that the nth moment has a scaling \(n(1-\gamma _a-\gamma _b)\), with intermittency correction \(\gamma _a(\frac{\beta _a^n-1}{\beta _a-1})+\gamma _b(\frac{\beta _b^n-1}{\beta _b-1})\). This is different from the three-dimensional stochastic Navier–Stokes equation (Birnir 2013a), where the moments are skewed Gaussians, but the structure functions have a scaling with intermittency corrections [we are making this comparison because this is the only other theory for a nonlinear stochastic PDE, with jumps, that has been worked out in similar details, see Birnir (2013a, 2013b).] The structure functions of Eq. (13) have the same scaling as the moments, however moments of Eq. (13) are even functions of \(\xi \) and reach their maximum at \(\xi = 0\), whereas the structure functions are odd functions of \(\xi _0\) and have two maxima at \(\xi _0 = \pm \text{ sech}^{-1}(\sqrt{3/2})\).
We thus see that the statistical theory of Eq. (13) 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
with \(B_0=1\). Thus, the (Brownian) noise in vessel branching only depends on the density (not all of its Fourier components) and is, in this aspect, similar to the noise in anastomosis, which also depends only on the jumps in the density. There is no reason to assume this at the beginning, so we assume that the noise is different for different (velocity) Fourier coefficients and obtain this rather surprising result.
Finally, the statistical theory of angiogenesis only applies to finite times, until a network of blood arteries or veins is beneficially established in organs and recovering tissue, or malignantly connected to 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 mimics the birth and death processes of tip branching and anastomosis, respectively, plus some other terms coming from the Brownian motion. Thus, our starting point is a stochastic PDE with Gaussian and Poisson noises that defines an infinite-dimensional Lévy process. We have solved the associated functional differential Kolmogorov–Hopf 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, in Sect. 5, the appropriate values of the parameters involved in the stochastic PDE.
Notes
We use boldface when we need to make clear that v and x are vectors, but otherwise indicate them in roman, with some abuse of notation.
References
Anderson, A.R.A., Chaplain, M.A.J.: Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60, 857–900 (1998)
Bauer, A.L., Jackson, T.L., Jiang, Y.: A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis. Biophys. J. 92, 3105–3121 (2007)
Bentley, K., Franco, C.A., Philippides, A., Blanco, R., Dierkes, M., Gebala, V., Stanchi, F., Jones, M., Aspalter, I.M., Cagna, G., Weström, S., Claesson-Welsh, L., Vestweber, D., Gerhardt, H.: The role of differential VE-cadherin dynamics in cell rearrangement during angiogenesis. Nat. Cell Biol. 16(4), 309–321 (2014)
Bernabeu, M.O., Jones, M.L., Nash, R.W., Pezzarossa, A., Coveney, P.V., Gerhardt, H., Franco, C.A.: PolNet: a tool to quantify network-level cell polarity and blood flow in vascular remodeling. Biophys. J. 114, 2052–2058 (2018)
Birnir, B.: The Kolmogorov–Obukhov statistical theory of turbulence. J. Nonlinear Sci. 23(4), 657–688 (2013)
Birnir, B.: The Kolmogorov–Obukhov Theory of Turbulence: A Mathematical Theory of Turbulence. Springer, New York (2013)
Bonilla, L.L., Capasso, V., Alvaro, M., Carretero, M.: Hybrid modeling of tumor-induced angiogenesis. Phys. Rev. E 90(6), 062716 (2014)
Bonilla, L.L., Carretero, M., Terragni, F.: Solitonlike attractor for blood vessel tip density in angiogenesis. Phys. Rev. E 94(6), 062415 (2016)
Bonilla, L.L., Carretero, M., Terragni, F., Birnir, B.: Soliton driven angiogenesis. Sci. Rep. 6, 31296 (2016)
Bonilla, L.L., Carpio, A., Carretero, M., Duro, G., Negreanu, M., Terragni, F.: A convergent numerical scheme for integrodifferential kinetic models of angiogenesis. J. Comput. Phys. 375, 1270–1294 (2018)
Bonilla, L.L., Carretero, M., Terragni, F.: Stochastic models of blood vessel growth. In: Giacomin, G., Olla, S., Saada, E., Spohn, H., Stoltz, G. (eds.) Stochastic Dynamics Out of Equilibrium. Springer Proceedings of Mathematics and Statistics, vol. 282, pp. 413–436. Springer, New York (2019)
Bonilla, L.L., Carretero, M., Terragni, F.: Two dimensional soliton in tumor induced angiogenesis. J. Stat. Mech. (2020). https://doi.org/10.1088/1742-5468/aba598
Byrne, H.M., Chaplain, M.A.J.: Mathematical models for tumour angiogenesis: Numerical simulations and nonlinear wave solutions. Bull. Math. Biol. 57, 461–486 (1995)
Byrne, H.M.: Dissecting cancer through mathematics: from the cell to the animal model. Nat. Rev. Cancer 10, 221–230 (2010)
Carmeliet, P.F.: Angiogenesis in life, disease and medicine. Nature 438, 932–936 (2005)
Carpio, Ana, Duro, Gema: Well posedness of an integrodifferential kinetic model of Fokker–Planck type for angiogenesis. Nonlinear Anal. Real World Appl. 30, 184–212 (2016)
Carpio, Ana, Duro, Gema, Negreanu, Mihaela: Constructing solutions for a kinetic model of angiogenesis in annular domains. Appl. Math. Model. 45, 303–322 (2017)
Chaplain, M.A.J., Stuart, A.: A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor. Math. Med. Biol. J. IMA 10, 149–168 (1993)
Chaplain, M.A.J.: The mathematical modelling of tumour angiogenesis and invasion. Acta. Biotheor. 43, 387–402 (1995)
Da Prato, Giuseppe: An Introduction to Infinite-Dimensional Analysis. Springer, Berlin (2006)
Da Prato, G., Zabczyk, J.: Ergodicity for Infinite Dimensional Systems. Lecture Notes of the London Mathematical Society, vol. 229. Cambridge University Press, Cambridge (1996)
Folkman, J.: Tumor angiogenesis: therapeutic implications. N. Engl. J. Med. 285, 1182–1186 (1971)
Folkman, J.: Tumor angiogenesis. Adv. Can. Res. 19, 331–358 (1974)
Franco, C.A., Jones, M.L., Bernabeu, M.O., Geudens, I., Mathivet, T., Rosa, A., Lopes, F.M., Lima, A.P., Ragab, A., Collins, R.T., Phng, L.-K., Coveney, P.V., Gerhardt, H.: Dynamic endothelial cell rearrangements drive developmental vessel regression. PLoS Biol. 13, e1002125 (2015)
Gebala, V., Collins, R., Geudens, I., Phng, L.-K., Gerhardt, H.: Blood flow drives lumen formation by inverse membrane blebbing during angiogenesis in vivo. Nat. Cell Biol. 18(4), 443–450 (2016)
Gerhardt, H., Golding, M., Fruttiger, M., Ruhrberg, C., Lundkvist, A., Abramsson, A., Jeltsch, M., Mitchell, C., Alitalo, K., Shima, D., Betsholtz, C.: VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. J. Cell Biol. 161(6), 1163–1177 (2003)
Heck, T., Vaeyens, M.M., Van Oosterwyck, H.: Computational models of sprouting angiogenesis and cell migration: towards multiscale mechanochemical models of angiogenesis. Math. Model. Nat. Phenomena 10, 108–141 (2015)
Hellström, M., Phng, L.K., Hofmann, J.J., Wallgard, E., Coultas, L., Lindblom, P., Alva, J., Nilsson, A.-K., Karlsson, L., Gaiano, N., Yoon, K., Rossant, J., Iruela-Arispe, M.L., Kalén, M., Gerhardt, H., Betsholtz, C.: Dll4 signalling through Notch1 regulates formation of tip cells during angiogenesis. Nature 445(7129), 776–780 (2007)
Hopf, Eberhard: Statistical hydromechanics and functional calculus. J. Ration. Mech. Anal. 1, 87–123 (1952)
Jackson, T., Zheng, X.: A cell-based model of endothelial cell migration, proliferation and maturation during corneal angiogenesis. Bull. Math. Biol. 72(4), 830–868 (2010)
Jolly, M.K., Boareto, M., Lu, M., Onuchic, J.N., Clementi, C., Ben-Jacob, E.: Operating principles of Notch–Delta–Jagged module of cell–cell communication. New J. Phys. 17, 055021 (2015)
Liotta, L.A., Saidel, G.M., Kleinerman, J.: Diffusion model of tumor vascularization. Bull. Math. Biol. 39, 117–128 (1977)
Mantzaris, N.V., Webb, S., Othmer, H.G.: Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 49, 111–187 (2004)
Martinson, W.D., Byrne, H.M., Maini, P.K.: Evaluating snail-trail frameworks for leader-follower behavior with agent-based modeling. Phys. Rev. E 102, 062417 (2020)
Martinson, W.D., Ninomiya, H., Byrne, H.M., Maini, P.K.: Comparative analysis of continuum angiogenesis models. J. Math. Biol. 82, 21 (2021)
Plank, M.J., Sleeman, B.D.: Lattice and non-lattice models of tumour angiogenesis. Bull. Math. Biol. 66, 1785–1819 (2004)
Øksendal, B.K., Sulem, A.: Applied Stochastic Control of Jump Diffusions, vol. 498. Springer, New York (2005)
Øksendal, B.: Stochastic differential equations. An Introduction with Applications. Springer, New York (2003)
Perfahl, H., Hughes, B.D., Alarcón, T., Maini, P.K., Lloyd, M.C., Reuss, M., Byrne, H.M.: 3D hybrid modelling of vascular network formation. J. Theor. Biol. 414, 254–268 (2017)
Pillay, S., Byrne, H.M., Maini, P.K.: Modeling angiogenesis: a discrete to continuum description. Phys. Rev. E 95, 012410 (2017)
Scianna, M., Bell, J., Preziosi, L.: A review of mathematical models for the formation of vascular networks. J. Theor. Biol. 333, 174–209 (2013)
Stokes, C.L., Lauffenburger, D.A.: Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J. Theor. Biol. 152, 377–403 (1991)
Stokes, C.L., Lauffenburger, D.A., Williams, S.K.: Migration of individual microvessel endothelial cells: stochastic model and parameter measurement. J. Cell Sci. 99, 419–430 (1991)
Szymborska, A., Gerhardt, H.: Hold me, but not too tight—endothelial cell–cell junctions in angiogenesis. Cold Spring Harb. Perspect. Biol. 10(8), a029223 (2018)
Terragni, F., Carretero, M., Capasso, V., Bonilla, L.L.: Stochastic model of tumour-induced angiogenesis: Ensemble averages and deterministic equations. Phys. Rev. E 93, 022413 (2016)
Travasso, R.D.M., Corvera Poiré, E., Castro, M., Rodríguez-Manzaneque, J.C., Hernández-Machado, A.: Tumor angiogenesis and vascular patterning: a mathematical model. PLoS ONE 6(5), e0019989 (2011)
Van Oers, R.F.M., Rens, E.G., La Valley, D.J., Reinhart-King, C.A., Merks, R.M.H.: Mechanical cell–matrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLoS Comput. Biol. 10(8), el003774 (2014)
Vega, R., Carretero, M., Travasso, R.D.M., Bonilla, L.L.: Notch signaling and taxis mechanism regulate early stage angiogenesis: a mathematical and computational model. PLoS Comput. Biol. 16(1), e1006919 (2020)
Vilanova, G., Colominas, I., Gomez, H.: Computational modeling of tumor-induced angiogenesis. Arch. Comput. Methods Eng. 24, 1071–1102 (2017)
Weinstein, Michael I.: Existence and dynamic stability of solitary wave solutions of equations arising in long wave propagation. Commun. Partial Differ. Equ. 12(10), 1133–1173 (1987)
Acknowledgements
This work was supported by a Chair of Excellence at the Universidad Carlos III de Madrid, in the spring of 2015, that is gratefully acknowledged by BB. LLB, MC and FT’s work has been supported by FEDER/Ministerio de Ciencia, Innovación y Universidades—Agencia Estatal de Investigación Grant No. PID2020-112796RB-C22, by the Madrid Government (Comunidad de Madrid, Spain) under the Multiannual Agreement with UC3M in the line of Excellence of University Professors (EPUC3M23) and in the context of the V PRICIT (Regional Programme of Research and Technological Innovation).
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Rustum Choksi.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix A: Jumps and Lévy Processes
In this appendix, we digress to explain the Ito–Lévy’s formula, a generalization of Ito’s formula that is necessary to solve Eq. (9). First we define stochastic processes with jumps, following Øksendal and Sulem (2005), where more information can be found. A Lévy process is a stochastic process on a filtered probability space \((\Omega ,{{\mathcal {F}}},\{{{\mathcal {F}}}_t\}_{t\ge 0}, {\mathbb {P}})\), which takes its values in \({\mathbb {R}}\), is continuous in probability and has stationary independent increments. The \(\{{{\mathcal {F}}}_t\}_{t\ge 0}\) is an increasing sequence of sigma-algebras indexed by t. They are called a filtration on the probability space \((\Omega , {{\mathcal {F}}}, {\mathbb {P}})\).
Let \(b_t\) be a one-dimensional Brownian motion on the probability space \((\Omega , {{\mathcal {F}}}, {\mathbb {P}})\) and suppose that N(t, z) is the number process of a Lévy process \(\eta _t\). An Ito–Lévy process is a stochastic process \(x_t\) on \((\Omega , {{\mathcal {F}}}, {\mathbb {P}})\) of the form
Here \({{\bar{N}}}\) is called the compensated jump measure of \(\eta _t\), defined as
or
where \(m(U) =\int _U E(N(\hbox {d}z,1))\) is the so-called Lévy measure of \(\eta _t\), see Birnir (2013b), Section 1.5.
The main computational tool in the Ito’s calculus for Ito–Lévy processes is the Ito–Lévy formula, see Øksendal and Sulem (2005). It is a generalization of the Ito’s formula, so that the jumps are included. Let \(x_t\) be the Ito–Lévy process
Let \(g(t,x)\in C^2([0,\infty ) \times {\mathbb {R}})\) be twice continuously differentiable. Then
is also an Ito–Lévy process and
Here \((udt+\sqrt{\sigma } db_t)^2\) has been computed by the rules
It is illustrative to use Eq. (A.2) to solve a differential equation. We now do so creating the geometric Lévy process. We solve the differential equation
Dividing by \(Z_t\)
we see that a reasonable guess for the function g above is
Applying Ito–Lévy formula (A.2), we get that
because \(\frac{\partial g}{\partial t} =\frac{\partial \ln (Z_t)}{\partial t}=0\). Integrating
with respect to t and exponentiating, we get that
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.
Appendix B: Moments from the Stochastic Model of Angiogenesis
In this appendix, we indicate how averages and moments are directly computed by numerical simulation of the stochastic angiogenesis model of Eq. (3) plus random branching and anastomosis processes. As in Terragni et al. (2016), we run \({\mathcal {N}}\) realizations of the underlying stochastic process and label each one by \(\omega \). Thus, \(\omega =1,\ldots ,{\mathcal {N}}\). The ‘microscopic marginal density’ is now
for realization \(\omega \), whereas the marginal density \(\tilde{p}(t,x)\) is the average of Eq. (B.1) over all realizations as \({\mathcal {N}}\rightarrow \infty \) and \(\sigma _x\rightarrow 0\). In practice, we select \({\mathcal {N}}\) so large that adding more realizations does not change the results of the simulations. Typically \({\mathcal {N}}=400\) suffices. Numerical values of all other parameters can be found in Terragni et al. (2016).
To calculate the nth moment of the marginal density, we use a similar formula
in which the ensemble average over realizations is explicitly written. Analogous calculations produce the structure function.
The numerical simulations of this paper consider the simple slab geometry of Terragni et al. (2016), Bonilla et al. (2016a), explained in Sect. 5. The soliton is calculated using data at a time when its development is completed to solve the collective coordinate equations for K and c in Eq. (22); see Bonilla et al. (2016a) for details.
Appendix C: Extension of the Mean-Value Theorem
We show that the integral
where p lies in the Hilbert space \(\mathcal {H}\), evaluates to \(\ln (p^*)\), the natural logarithm of a positive constant \(p^*\), by a slight extension of the mean-value theorem for integrals.
First, we recall from the end of Sect. 4 that \(p_\infty =p_\infty (\xi )\) is a traveling wave given by the soliton. Thus, all the directions except that given by the soliton integrate to one in the above integral and we get that
where \(e_s\) and \(q_s\) are the mean and variance associated with the direction of the soliton in the function space. However, as \(\xi \) ranges from \(-\infty \) to \(\infty \), p ranges from 0 to its maximum, and from its maximum to 0. By the chain rule, this amounts to integrating the maximum, \(p=\hbox {max}\ p(\xi )= \max \ p_\infty (\xi )\), twice from 0 to \(\infty \), namely
If the interval was finite, i.e., [0, C], we immediately get that
where \(p^*\) is a positive constant, by the mean-value theorem for integrals. We use the following argument to reach the same conclusion for \(C=\infty \). First, the integral in the left-hand side is clearly increasing as a function of C, for C large. Then, since
\( \ln (p^*(C))\) is increasing as well. However, \(\ln (p^*(C))\) is also bounded as
and
for C large enough and some \(\epsilon > 0\), so
Now, a bounded increasing sequence of numbers must converge and we set
Then
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Birnir, B., Bonilla, L., Carretero, M. et al. The Statistical Theory of the Angiogenesis Equations. J Nonlinear Sci 34, 29 (2024). https://doi.org/10.1007/s00332-023-10006-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00332-023-10006-2
Keywords
- Angiogenesis
- Stochastic PDEs
- Statistical theory of nonlinear PDEs
- Moments and structure functions
- Mathematical biology