A coupled stochastic differential reaction-diffusion system for angiogenesis

A coupled system of nonlinear mixed-type equations modeling early stages of angiogenesis is analyzed in a bounded domain. The system consists of stochastic differential equations describing the movement of the positions of the tip and stalk endothelial cells, due to chemotaxis, durotaxis, and random motion; ordinary differential equations for the volume fractions of the extracellular fluid, basement membrane, and fibrin matrix; and reaction-diffusion equations for the concentrations of several proteins involved in the angiogenesis process. The drift terms of the stochastic differential equations involve the gradients of the volume fractions and the concentrations, and the diffusivities in the reaction-diffusion equations depend nonlocally on the volume fractions, making the system highly nonlinear. The existence of a unique solution to this system is proved by using fixed-point arguments and H\"older regularity theory. Numerical experiments in two space dimensions illustrate the onset of formation of vessels.


Introduction
Angiogenesis is the process of expanding existing blood vessel networks by sprouting and branching. Its mathematical modeling is important to understand, for instance, wound healing, inflammation, and tumor growth. In this paper, we analyze a variant of the continuum-scale model suggested in [2] that is used to simulate the early stages of angiogenesis. The model takes into account the dynamics of the tip (leading) endothelial cells by solving stochastic differential equations, the influence of various proteins triggering the cell dynamics by solving reaction-diffusion equations, and the change of some components of the extracellular matrix into extracellular fluid by solving ordinary differential equations. Up to our knowledge, this is the first analysis of the model of [2].
Angiogenesis is mainly triggered by local tissue hypoxia (low oxygen level in the tissue), which activates the production of the signal protein vascular endothelial growth factor (VEGF). Endothelial cells, which form a barrier between vessels and tissues, reached by the VEGF signal initiate the angiogenic program. These cells break out of the vessel wall, degrade the basement membrane (a thin sheet-like structure providing cell and tissue support), proliferate, and invade the surrounding tissue while still connected with the vessel network. The angiogenic program specifies the activated endothelial cells into tip cells (cells at the front of the vascular sprouts) and stalk cells (highly proliferating cells). The tip cells lead the sprout towards the source of VEGF, while the stalk cells proliferate to follow the tip cells supporting sprout elongation; see Figure 1. Following [2], the tip cells segrete the proteins delta-like ligand 4 (DLL4), matrix metalloproteinase (MMP), and urokinase plasminogen activator (uPA). The chemokine DLL4 makes the stalk cells follow the tip cells, MMP breaks down the basement membrane, and uPA degrades the fibrin matrix such that cells can move into. There are many other molecular mechanisms and mediators in the angiogenesis process; see, e.g., [1,28] for details. where t ≥ 0 is the time and x ∈ D ⊂ R 3 the spatial variable. All unknowns depend additionally on the stochastic variable ω ∈ Ω, where Ω is the set of events. We assume that the mixture of basement membrane, extracellular fluid, and fibrine matrix is saturated, i.e., the volume fractions f B , f E , and f F sum up to one. We introduce the vectors X = (X k i ) i=1,2,k=1,...,N i , f = (f B , f E , f F ), and c = (c V , c D , c M , c U ).
Stochastic differential equations. The tip cells move according to chemotaxis force, driven by the gradient of the VEGF concentration, the durotaxis force, driven by the gradient of the solid fraction f S := f B +f F , and random motion modeling uncertainties. The dynamics of X k i (t) is assumed to be governed by the stochastic differential equations (SDEs) are Wiener processes, and the drift terms (2) g 1 [c, f ](X 1 , t) := α 0 M 1 (f S (X 1 ), X 1 )z 1 + γ(f S (X 1 ))∇c V (X 1 , t) + λ(f S (X 1 ))∇f S (X 1 , t), include a constant α 0 > 0, the strain energy M i , the direction of movement determined by the strain energy density z i , the chemotaxis force ∇c V in the direction of VEGF (tip cells) and ∇c D in the direction of DLL4 (stalk cells), and the migration as a result of the durotaxis force ∇f S . In this paper, we allow for general drift terms by imposing suitable Lipschitz continuity conditions; see Assumption (A4) below. In the numerical experiments, we choose the functions M i , γ, and λ as in Appendix B.
Ordinary differential equations. The proteins MMP and uPA degrade the basement membrane and fibrin matrix, respectively, while enhancing the extracellular fluid component. Therefore, the volume fractions f B , f E , and f F are determined by the ordinary differential equations (ODEs) where s B , s F > 0 are some rate constants. Note that the last differential equation is redundant because of the volume-filling condition f B + f E + f F = 1. Clearly, equations (3) can be solved explicitly, giving for (x, t) ∈ D × (0, T ) and pathwise in Ω (we omit the argument ω), Reaction-diffusion equations. The mass concentrations are modeled by reaction-diffusion equations, describing the consumption and production of the proteins: with initial and no-flux boundary conditions where the rate terms are given by for j = V, D, M, U and V k j : Ω × R 3 → R are nonnegative smooth potentials approximating the delta distribution. The parameters r j and s j are positive. In [2], the rate terms are given by delta distributions instead of smooth potentials. We assume smooth potentials because of regularity issues, but they can be given by delta-like functions as long as they are smooth. Indeed, we need C 1+δ (D) solutions c j to solve the SDEs (1), and this regularity is not possible when the source terms of (5) include delta distributions. As the number of the proteins is typically much larger than the number of tip cells, the stochastic fluctuations in the concentrations are expected to be much smaller than those associated with the tip cells, which justifies the macroscopic approach using reaction-diffusion equations for the concentrations.
In equation (5) for c V , the term α V c V models the consumption of VEGF along the trajectory of the tip cells. The protein DLL4 is regenerated from conversion of VEGF, modeled by α D c V along the trajectories of the tip cells, and consumed by the stalk cells, modeled by β D c D along the trajectories of the stalk cells. In equation (5) for c M , the term s M f B c M describes the decay of the MMP concentration with rate s B > 0 as a result of the breakdown of the basement membrane, and α M c V models the production of MMP due to conversion from VEGF. Similary, s U f F c U describes the decay of the uPA concentration, which breaks down the fibrin matrix, and the protein uPA is regenerated, leading to the term α U c V .
The diffusivities are given by the mixing rule where D i j > 0 for i = B, E, F . Then and equations (5) are uniformly parabolic. Note, however, that equations (5) are nonlocal and quasilinear, since the diffusivities are determined by the time integrals of c M or c U ; see (4). Various biological phenomena are not modeled by our equations. For instance, we do not include the initiation of sprouting from preexisting parental vessels, the branching from a tip cell, and anastomosis (interconnection between blood vessels). Moreover, in contrast to [2], we do not allow for the transition between the phenotypes "tip cell" and "stalk cell" to simplify the presentation. On the other hand, we may include further angiogenesis-related proteins, if the associated reaction-diffusion equations are of the structure (5).
1.2. State of the art. There are several approaches in the literature to model angiogenesis, mostly in the context of tumor growth. Cellular automata models divide the computational domain into a discrete set of lattice points, and endothelial cells move in a discrete way. Such models are quite flexible, and intra-cellular adherence can be easily implemented, but their numerical solution is computationally expensive when the numbers of cells or molecules are large [12]. In semicontinuous models, the cells are treated as discrete entities, but their movement is not restricted to any lattice points [4]. Continuum-scale models consider cells as averaged quantities, leading to cell densities whose dynamics is described by partial differential equations; see, e.g., [3] for wound healing and [10] for angiogenesis. Chemotaxis can be modeled in this approach by Keller-Segel-type equations, which admit global weak solutions in two space dimensions without blowup [9]. A hybrid approach was investigated in [5], where the blood vessel network is implemented on a lattice, tip cells are moving in a lattice-free way, and other cells are modeled macroscopically as densities. A continuum-scale approach was chosen in [2], which is the basis of the present paper. The novelty of [2] is the distinction of tip and stalk cells and the inclusion of proteins secregated by the tip cells.
In other models, stochastic effects have been included. In [26], the movement of the tip cells is modeled by a SDE, with a deterministic part describing chemotaxis, and a stochastic part modeling random motion. The mean-field limit in a stochastic many-particle system, leading to reaction-diffusion equations, was performed in [6,25]. We also refer to the reviews [8,24] on further modeling approaches of angiogenesis.
Numerical simulations of a coupled SDE-PDE model for the movement of the tip cells and the dynamics of the tumor angiogenesis factor, fibronectin (a protein of the extracellular matrix), and matrix degrading enzymes were presented in [7]. Other SDE-PDE models in the literature are concerned with the proton dynamics in a tumor [17], acid-mediated tumor invasion [13], and viscoelastic fluids [14]. However, only the works [13,17] treat a genuine coupling between SDEs and PDEs. While the model in [13] also includes a cross-diffusion term in the equation for the cancer cells, we have simpler reaction-diffusion equations but with nonlocal diffusivities. Up to our knowledge, the mathematical analysis of system (1), (3)-(6) is new.
1.3. Assumptions and main result. Let (Ω, F, (F t ) t≥0 , P) be a stochastic basis with a complete and right-continuous filtration, and let (W k i (t)) t≥0 for i = 1, 2, k = 1, . . . , N i be independent standard Wiener processes on R d (d ∈ N) relative to (F t ) t≥0 . We write L p (Ω, F; B) for the set of all F-measurable random variables with values in a Banach space B, for which the L p norm is finite. Furthermore, let D ⊂ R 3 be a bounded domain with boundary ∂D ∈ C 3 (needed to obtain parabolic regularity; see Theorem 22). We set Q T = D × (0, T ). We write C k+δ (D) with k ∈ N 0 , δ ∈ (0, 1) for the space of C k functions u such that the kth derivative D k u is Hölder continuous of index δ. The space of Lipschitz continuous functions on D is denoted by C 0,1 (D). For notational convenience, we do not distinguigh between the spaces C k+δ (D; R n ) and C k+δ (D). Furthermore, we usually drop the dependence on the variable ω ∈ Ω in the ODEs and PDEs, which hold pathwise P-a.s. Accordingly, we write c instead of c(ω, ·, ·) and c(t) instead of c(ω, ·, t). Finally, we write "a.s." instead of "P-a.s.".
We impose the following assumptions: (A1) Initial data: at most linear growth in x, is progressively measurable, and satisfies σ i = 0 on ∂D a.s.
Furthermore, for c, f ∈ L ∞ (0, T ; W 2,∞ (D)), there exists L 2 > 0 such that for (x, t), (x , t) ∈ D × [0, T ] and i = 1, 2, Let us discuss these assumptions. Since we need C 1+δ solutions (c, f ) to obtain Hölder continuous coefficients of the SDEs (which ensures their solvability), we need some regularity conditions on the initial data in Assumption (A1). Accordingly, ∇c 0 j · ν = 0 on ∂D is a compatibility condition needed for such a regularity result. In Assumption (A1), we impose the volume-filling condition initially, f 0  (1). Note that g i [c, f ] in example (2) satisfies Assumption (A4). Finally, Assumption (A5) is a simplification to ensure the parabolic regularity results needed, in turn, for the solvability of (1).  (21) below. Then α j and β j are Hölder continuous in D × [0, T ] a.s. as a function of X.
As the first step, we prove the solvability of the ODEs (3) with c j only lying in the space L ∞ (0, T ; L 2 (D)). This is achieved by an approximation argument, assuming a sequence c (k) j ∈ L ∞ (Q T ) and then passing to the limit. This is possible since the solutions f i , represented in (4), are bounded.
Second, we show that there exists a weak solution c, which is defined pathwise for ω ∈ Ω, to the reaction-diffusion equations (5) by using a Galerkin method and the fixed-point theorem of Schauder for fixed α j and β j . The regularity of c can be improved step by step.
Moser's iteration method shows that the weak solution is in fact bounded, and a general regularity result for evolution equations yields ∂ t c ∈ L 2 (Q T ). Then, interpreting equations (5) as elliptic equations with right-hand side ∂ t c ∈ L 2 (Q T ), we conclude the Hölder continuity of c(t) for any fixed t ∈ (0, T ). Thus, the diffusivities are Hölder continuous, and we infer C 1+δ (D) and W 2,∞ (D) solutions c. We conclude classical solvability from the regularity results of Ladyženskaya et al. [18].
In the third step, we solve the SDEs (1). The functions (c, f ) have Hölder continuous gradients, and we show that (c, f ) and (∇c, ∇f ) are progressively measurable. Therefore, together with Assumption (A4), the conditions of the existence theorem of [21, Theorem 3.1.1] are satisfied, and we obtain a solution X to (1) in the sense (9).
Fourth, defining a suitable complete metric space Y R (0, T ; D), endowed with the C 0 ([0, T ]; L 4 (D)) norm, this defines the fixed-point operator Φ : X → X on Y R (0, T ; D). The fixedpoint operator can be written as the contatenation . It remains to show that Φ is a self-mapping and a contraction, which is possible for a sufficiently large R > 0. In fact, we show that for any n ∈ N, for all X, X ∈ Y R (0, T ; D), where c n → 0 as n → ∞. It follows from the variant of Banach's fixed-point theorem in [19,Theorem 2.4] that Φ has a unique fixed point, which gives a unique solution (X, c, f ) to (1), (3)- (7).
The paper is organized as follows. The solvability of the ODEs (3) for L 2 (D) coefficients and a stability estimate are shown in Section 2. In Section 3, the existence of a unique classical solution to the reaction-diffusion equations (5)-(6) and some stability results are proved. The progressive measurability of the solutions to (3) and (5) as well as the solvability of the SDEs (1) is verified in Section 4. Based on these results, Theorem 1 is proved in Section 5. Some numerical experiments are illustrated in Section 6, showing stalk cells following the tip cells and forming a premature sprout. Appendix A summarizes some regularity results for elliptic and parabolic equations used in this work, and Appendix B collects the numerical values of the various parameters used in the numerical experiments.

Solution of ODEs in Bochner spaces
We prove a result for ODEs in some Bochner space, which is needed for the solution of (3) when the concentrations c j are only L 2 (D) functions.
. This proves the claim. Next, we consider the differential equation . We want to prove that (g k ) converges as k → ∞ to a solution of the original equation. It follows that We conclude from Gronwall's lemma and (11) that Thus, (g k (t)) is a Cauchy sequence for every t ∈ [0, T ] and we infer that g k (t) → g(t) in L 2 (D). Because of the L ∞ (0, T ; L 2 (D)) bound for (g k ) and dominated convergence again, ). There exists a subsequence, which is not relabeled, such that g k (t) → g(t) a.e. in D, for any t ∈ [0, T ]. We deduce from g k (t) ≤ K that g(t) ≤ K for all t ∈ [0, T ]. This shows that g ∈ L ∞ (Q T ). Writing (12) as an integral equation and performing the limit k → ∞, the previous convergence results show that g solves (10).
We also need the following stability result, which relates the difference f 1 −f 2 of solutions to (3) with the difference c 1 − c 2 of solutions to (5).
Then there exists C > 0, only depending on T , the L ∞ (0, T ; W k,∞ (D)) norm of u i , and the W k,∞ (D) norm of g 0 i , such that for p > 1, ) . Proof. The regularity of g i follows from the explicit formula and the regularity for g 0 i and u i . Furthermore, taking the W k,p (D) norm of the result follows from the regularity u 1 ∈ L ∞ (0, T ; W k,∞ (D)) and g 2 ∈ W k,∞ (D).

Solution of the reaction-diffusion equations
We show the existence of a weak solution to (5)-(6) for given α, β D and then prove that this solution is in fact a classical one. Lemma 4 (Existence for (5) Proof. Let (e ) ∈N be the sequence of eigenfunctions of the Laplace operator with homogeneous Neumann boundary conditions. They form an orthonormal basis of L 2 (D) and are orthogonal with respect to the H 1 (D) norm. The boundary regularity (13), we consider first a linearized version. Let c N ∈ C 0 ([0, T ]; E N ) be given. Then  (15) c N L ∞ (0,T ;L 2 (D)) + c N L 2 (0,T ;H 1 (D)) + ∂ t c N L 2 (0,T ;H 1 (D) ) ≤ C(T, c 0 ), c N L ∞ (0,T ;H 1 (D)) ≤ C(N ). In the following, the notation (E N ; · Z ) means that E N is endowed with the norm of the Banach space Z, requiring that E N ⊂ Z.
We wish to define a fixed-point operator c N → c N . For this, we introduce the sets and W is defined as the closure of Y 4 b with respect to the C 0 ([0, T ]; L 2 (D)) norm. We claim that the embedding W → C 0 ([0, T ]; (E N ; L 2 (D))) is compact. Indeed, (E N ; H 1 (D)) is compactly embedded into (E N ; L 2 (D)). By the Aubin-Lions lemma, Y is compactly This defines the fixed-point operator S : W → C 0 ([0, T ]; (E N ; L 2 (D))). Estimates (15) show that S(W ) ⊂ Y 4 b ⊂ W . Furthermore, we can prove by standard arguments that S : W → W is continuous. We deduce from Schauder's fixed-point theorem that there exists a fixed point c N of S satisfying (13). Now, we perform the limit N → ∞. Estimates (15) allow us to apply the Aubin-Lions lemma, giving the existence of a subsequence (not relabeled) such that, as N → ∞, .
We define We infer from the linear dependence of . Thus, passing to the limit N → ∞ in equation (13) and similarly for the other equations in (13). Here, ·, · denotes the dual product between H 1 (D) and H 1 (D) . Since the class of functions φ of this type is dense in L 2 (0, T ; H 1 (D)) ∩ H 1 (0, T ; H 1 (D) ) [30,Prop. 23.23iii], we conclude that (c, f ) is a weak solution to (3)- (6).
It remains to verify that the limit concentration c is nonnegative componentwise. We use c − V = min{0, c V } as a test function in the weak formulation of equation (5) for c V : which yields c − V (t) = 0 and c V (t) ≥ 0 in D, t > 0. Next, we use c − D as a test function in the weak formulation of equation (5) for c D : since c V ≥ 0 and c − D ≤ 0. We infer that c D (t) ≥ 0 in D, t > 0. In a similar way, we prove that c M (t) ≥ 0 and c U (t) ≥ 0 in D, t > 0. This finishes the proof.

3.2.
Regularity. The existence theory for SDEs requires more regularity for the solution c to (5). First, we prove L ∞ bounds. We suppose that Assumptions (A1)-(A5) hold in this section.
Lemma 5. Let c be a weak solution to (5)- (6). Then c ∈ L ∞ (Q T ) and it holds for all as a test function in the weak formulation of equation (5) for c V : . Then (c D (0) − M ) + = 0 and, choosing (c D − M ) + as a test function in the weak formulation of equation (5) for c D , where we used β D ≥ 0, and the last inequality follows from the choice of M 0 . This shows that c D (t) ≤ M 0 e t in D, t > 0. The bounds for c M and c U are shown in an analogous way.
Next, we prove that the solution c(t) is Hölder continuous.
Lemma 6. Let c be a weak solution to (5)- (6). We suppose that there exists Λ > 0 such that for a.e. 0 < t < T , Then there exists δ > 0 such that for 0 < t < T , where C 2 > 0 depends on the L 2 (D) norm of c 0 , the L ∞ (D) norm of f 0 , and Λ, and δ, C δ depend on the lower and upper bounds (8) for D j and Λ.
Proof. The L 2 (Q T ) bound for ∂ t c follows immediately from Theorem 19 in the Appendix. The Hölder estimate follows from [23,Prop. 3.6]. Indeed, we interpret equation (5) for as an elliptic equation with bounded diffusion coefficient D V (f ) and right-hand side in L p (D) with p > d/2. By [23,Prop. 3.6], there exists δ > 0 such that c V (t) ∈ C δ (D) and The result follows by observing that d ≤ 3 implies that p < 2. The dependency of δ and C δ on the data follows from [11,Theorem 8.24], which is the essential result needed in the proof of [23,Prop. 3.6]. The regularity for the other concentrations is proved in a similar way.
Lemma 7. Let c be a weak solution to (5)- (6). Furthermore, let c 0 Proof. We know from Lemma 6 that c(t) is Hölder continuous in D for a.e. t ∈ (0, T ). We claim that f is Hölder continuous in D × [0, T ]. Let x, y ∈ D and τ, t ∈ [0, T ]. We assume without loss of generality that τ < t. The Lipschitz continuity of z → exp(−z) implies, using the explicit formula for f B , that , where we also used Lemma 6. Similar estimates hold for f E and f F . Thus, the assumptions of Theorem 20 in the Appendix are fulfilled, yielding the statement.
For the solvability of the SDEs (1), we need the W 2,∞ (D) regularity of c. For this, we introduce the following space: where q ∈ [1, ∞], equipped with the norm Proof. We show the statement for c V , as the proofs for the other components are similar. f (x, t))∇c 0 (x)) is bounded by a constant for (x, t) ∈ D × (0, T ). This constant depends on Λ, c 0 , f 0 , and T but not on c. By Lemma 7, the diffusion coefficient D V (f ) is an element of C 1+δ,(1+δ)/2 (D × [0, T ]). Therefore, an application of Theorem 21 finishes the proof.
For smooth initial data, we can obtain classical solutions to (5)- (6). The following result is not needed for the solvability of (1) but stated for completeness.
Proof. The result follows from Theorem 22 in the Appendix after bringing (5) in nondivergence form.
3.3. Stability and uniqueness. The stability results are used for the solution of the SDEs; they also imply the uniqueness of weak solutions. We start with a stability estimate in the norms L ∞ (0, T ; L 2 (D)) and L 2 (0, T ; H 1 (D)). Let Assumptions (A1)-(A5) hold.
Lemma 10. Let c i for i = 1, 2 be weak solutions to (5)-(6) with the same initial data (c 0 , f 0 ) but possibly different coefficients α i and β i . Then there exists C > 0, which is independent of c i , such that for all t ∈ [0, T ], Proof. We first consider c V . We take the difference of the equations satisfied by c 1,V − c 2,V and take the test function c 1,V − c 2,V in its weak formulation. This leads to 1 2 Let ε := min{D i j : j = V, D, M, U, i = B, E, F } > 0. Using Young's inequality and the estimate (f 1 − f 2 )(t) L 2 (D) ≤ C c 1 − c 2 L 1 (0,T ;L 2 (D)) from Lemma 3, we find that . The estimates for c 1,j − c 2,j with j = D, M, V are similar. This gives d dt . An application of Gronwall's lemma finishes the proof.
A stability estimate can also be proved with respect to the H 2 (D) norm.
Proof. The difference u := c 1,V − c 2,V is the solution to the linear problem where, by Lemma 8, the right-hand side is an element of L 2 (Q T ). Since the diffusion coefficient is bounded, we can apply Theorem 19 in the Appendix to conclude that For the estimate of the right-hand side, we recall from Lemma 3 that Then, using the linearity of D V and the estimate for c 2,V from Lemma 8, we infer that The difference c 1 − c 2 in the L 2 (0, T ; H 1 (D)) norm can be estimated according to Lemma 10. Therefore, Similar estimates can be derived for the differences c 1,j − c 2,j (j = D, M, U ).
To estimate u in the L 2 (0, T ; H 2 (D)) norm, we use the inequality Thus, it remains to consider ∆u. We deduce from We infer from (19) and related inequalities for c 1,j − c 2,j that which concludes the proof.
Proof. Let u = c 1,V − c 2,V be the solution to (17). Since g ∈ L 4 (Q T ) by Lemma 8 (recall definition (18) of g), Theorem 21 in the Appendix shows that u ∈ W 2,1,4 (Q T ) and, because of ∇c 2,V ∈ L ∞ (Q T ), The first term on the right-hand side can be estimated by using the embedding H 2 (D) → W 1,4 (D) and Lemma 11: We deduce from the embedding H 1 (D) → L 4 (D) and Lemma 11 again that This gives The estimates for c 1,j − c 2,j (j = D, M, U ) are similar.

Solution of the stochastic differential equations
Let α, β be given by (7). We first study the measurability of (c, f ).
Lemma 13. Let f 0 ∈ L ∞ (Ω; C 1+δ (D)) and c 0 ∈ L ∞ (Ω; W 2,∞ (D)) such that ∇c 0 j · ν = 0 on ∂D, j = V, D, M, U . Furthermore, let (c, f ) be a pathwise solution to Proof. Since f j can be represented as a function depending on the time integral of c, it is sufficient to show the measurability of c j . The continuity of the potentials defining α j and β j in (7) shows that α j and β j are processes with càdlàg paths almost surely. By approximating the initial data c 0 , f 0 and the processes α j , β j by suitable simple processes, which are adapted to the filtration by construction, we can obtain the F t -measurability of c j (t) : Ω → C 1 (D) for t ∈ [0, T ]. We conclude from Lemma 10 and the compactness of W 2,∞ (D) ⊂ C 1+δ (D) in C 1 (D) the measurability of c j as the limit of measurable functions. For details of this construction, we refer to [13,Section 3.3]. It is known that càdlàg processes, which are adapted to the filtration, are also progressively measurable [16,Prop. 1.13]. If the filtration is complete, this holds also true for processes having càdlàg paths almost surely. The estimate c j (t) − c j (s) C 1 (D) ≤ C|t − s| δ , which follows from Lemma 6, implies that c j (t) has almost surely continuous paths and consequently, c j (t) is progressively measurable. To be precise, this yields the measurability of c j as a function from (Ω × [0, t], F t × B([0, t])) to (C 1 (D), B(C 1 (D))) for every t ∈ [0, T ].
The function (c, x) → c(x), C 1 (D) × D → R 4 , is continuous and hence, it is measurable as a mapping from (C 1 (D), B(C 1 (D))) to (R 4 , B(R 4 )). Now, we can write c(ω, x, t) as the concatenation of measurable functions, which yields the measurability of c. In a similar way, we can prove the measurability of ∂c j /∂x i for i = 1, 2, 3 by considering the continuous mapping (c, x) → (∂c/∂x i )(x).

Lemma 14.
Let Assumptions (A1)-(A5) hold. Then there exists a unique, progressively measurable solution (X k i ) to (1) such that X k i (t) ∈ D a.s. for every t ∈ [0, T ]. Proof. We extend the coefficients g i and σ i by setting them to zero outside of D. The extended coefficients are still uniformly Lipschitz continuous. We infer from Lemma 13 that g i is progressively measurable. Thus, by [15,Theorem 32.3], there exists a strong solution to (1).
It remains to show that X k i (t) ∈ D a.s. Let φ be a smooth test function satisfying supp φ ⊂ D c . We obtain from Itô's lemma that (20) , t) = 0 by Assumption (A3) and σ i (X k i (t)) = 0 by Assumption (A2). Equation (20) then shows that φ(X k i (t)) = φ(X 0 i ) = 0 and X k i (t) ∈ (supp φ) c a.s. Since φ with supp φ ⊂ D c was arbitrary, we conclude that X k i (t) ∈ D a.s. for t ∈ (0, T ).

Proof of Theorem 1
The fixed-point operator is defined as a function that maps X → (α, β D ) → (c, f ) → X, where (α, β D ) are defined in (7) with X replaced by X. To define its domain, we need some preparations. For given R > 0, we introduce the following space: Proof. Let (X n ) be a Cauchy sequence in Y R (0, T ; D) and let ε > 0. Then there exists N ∈ N such that for all n, m ≥ N , For any t ∈ [0, T ], (X n (t)) is a Cauchy sequence in L 4 (Ω). Consequently, X n (t) → X(t) in L 4 (Ω), where X(t) ∈ L 4 (Ω) is F t -measurable. Furthermore, there exists a subsequence of (X n (t)) (not relabeled) that converges pointwise to X(t) a.s., proving that X(t) ∈ D a.s. The definition of the Hölder norm implies that X n (t) − X n (s) L 4 (Ω) ≤ R|t − s| 1/2 for all s, t ∈ [0, T ]. This gives in the limit n → ∞ that X(t) − X(s) L 4 (Ω) ≤ R|t − s| 1/2 and consequently X ∈ C 1/2 (0, T ; D). We conclude that X ∈ Y R (0, T ; D). By the Kolmogorov continuity criterium, (a modification of) X has almost surely Hölder continuous paths. As X(t) is an adapted process with respect to the filtration F t , X is progressively measurable.
Lemma 16. Let X ∈ Y R (0, T ; D) for some R > 0, and let (c, f ) be a solution to (3)-(6), where α, β are given by (7) with X replaced by X. Then there exists R 0 > 0 not depending on R such that the solution X to (1) satisfies X ∈ Y R 0 (0, T ; D).
Proof. According to Lemma 8, c is bounded in the L ∞ (0, T ; W 2,∞ (D)) norm by a constant that is independent of R. Then, by Lemma 14, there exists a unique solution X to (3). Since X k i (t) ∈ D a.s., c, ∇c, f , ∇f are bounded uniformly in R, i.e., there exists K = K(c 0 , f 0 ) > 0, which is independent of R, such that |g i [c, f ](X k i (t), t)| ≤ K a.s. Thus, for s, t ∈ [0, T ], using the Burkholder-Davis-Gundy inequality, The lemma follows after choosing R 0 := C(K, T, σ, D) 1/4 .
The previous lemma shows that the fixed-point operator Φ : Y R 0 (0, T ; D) → Y R 0 (0, T ; D), X → X, is well defined. We need to verify that Φ is a contraction. We first prove an auxiliary result.  6), where α, β are given by (7) with X replaced by X, X ∈ Y R (0, T ; D) for some R > 0, respectively. Then the associated solutions X and X to (1) satisfy where the constant C > 0 does not depend on R, (c, f ), or (c , f ).
Proof. The Itô integral representation of X(t) − X (t) gives It follows from Assumption (A4) that Furthermore, by the Burkholder-Davis-Gundy inequality and the Lipschitz continuity of σ, We insert these estimates into (22) and use Hölder's inequality: Then Gronwall's lemma concludes the proof.

Numerical experiments
We illustrate the dynamics of the tip and stalk cells in the two-dimensional ball D = B R (0) around the origin with radius R = 500 (in units of µm). Let h = 10 be the space step size and introduce the grid points x ij = ((k −i)h, (k −j)h) ∈ R 2 , where i, j = 0, . . . , 2k and k = R/h. The time step size equals τ = 1 (in units of seconds).
The stochastic differential equations (1) are discretized by using the Euler-Maruyama scheme. The nonlinearity g i [c, f ] is chosen as in (2) with M , γ, and λ given in Appendix B. Furthermore, α 0 and z are taken as in [29, formulas (10) and (14)]. Compared to [2], we neglect the contribution of the Hertz contact mechanics regarding z to guarantee the boundary condition g i [c, f ](·, t) = 0 on ∂D. We choose the continuous radially symmetric stochastic diffusion for 9R/10 < |x| < R, 1/10 for |x| ≤ 9R/10.
The solutions (4) to the ordinary differential equations (3) are written iteratively as and similarly for f F . The integral is approximated by the trapezoid rule where c n M,ij approximates c M (x ij , nτ ). We set f n ij := (f B , f E , f F )(x ij , nτ ). Finally, we discretize the reaction-diffusion equations (5) using the forward Euler method and the central finite-difference scheme . Notice that we obtain a semi-implicit scheme. The resulting linear system of equations is implemented in the Python-based software environment SciPy using sparse matrices and solved by using the spsolve function from the scipy.sparse.linalg package.
The potentials V k j , used in (5), are given by x ∈ D, j = D, M, U, V, where R m = 12.5, and I > 0 is a normalization constant to ensure that R 2 V k j (x)dx = 1. It remains to define the initial conditions. The initial positions of the endothelial cells X 0,k i (i = 1, 2, k = 1, . . . , N i ) are given by where (r, φ) is uniformly drawn from the set [0.65R, 0.75R] × [0, π/2] and R0.65 = 325, R c = 0.75R = 375. The initial volume fractions are as well as f 0 which is concentrated at the origin, and assume that the concentrations of the remaining proteins vanish, c 0 D = c 0 M = c 0 U = 0 in D, as they are segregated by the tip cells.
We choose N 1 = 2 tip cells and N 2 = 200 stalk cells. Figure 2 shows the positions of the tip and stalk cells at different times. The tip cells segregate the DLL4 protein, and the stalk cells detect the local increase of the DLL4 concentration, such that they follow the corresponding tip cell. This effect is slightly more pronounced for the tip cell that starts in an environment with a dense stalk cell population. The position of this tip cell is closer to the origin than the other tip cell with a higher VEGF concentration, leading to a relatively high production of DLL4 proteins. The stalk cells, which do not follow a tip cell, are primarily influenced by the stiffness gradient ∇(f B + f F ) and the strain energy density M , which incorporates contact mechanics, resulting to a spreading of these cells. The protein concentrations are shown in Figure 4. As the diffusion coefficient for VEGF is much larger than the reaction rate s V , the concentration of the VEGF protein becomes uniform in the large-time limit. The DLL4, MMP, and uPA proteins are produced by the tip cells and hence follow their paths. The corresponding concentrations increase with the availability of VEGF and decrease due to consumption by the stalk cells or by getting exhausted from breaking down the fibrin matrix or the boundary membrane. Since the diffusion is slow, the changes in the concentration are local up to time T = 1600 s.
We present the volume fractions of the basement membrane, fibrin matrix, and extracellular fluid in Figure 3. The membrane and fibrin matrix are degraded by the MMP and uPA proteins, thus increasing the volume fraction of the extracellular fluid. As both proteins are produced by the tip cells, the degradation follows their paths.
Summarizing, we see that the model successfully describes the formation of premature sprouts. The experiments from [2] for dermal endothelial cells show that the in vitro angiogenesis sprouting qualitatively well agrees with the numerical tests. Clearly, the proposed system of equations models only a very small number of biological processes, chemical reactions, and signal proteins, and more realistic results can be only expected after taking into account more biological modeling details. Still, the onset of vessel formation is well illustrated by our simple model.   Then u ∈ H k+2 (D), and there exists a constant C > 0 not depending on u or f such that u H k+2 (D) ≤ C f H k (D) + u L 2 (D) .

Appendix B. Model parameters and constants
The model parameters and constants are taken from [2]. For the convenience of the reader, we collect here the expressions: The parameters are chosen as in the following table; see [2,Appendix].